diff --git a/ccpp/physics b/ccpp/physics index 129183136..4fe67293f 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 129183136d20db972138609a23c0286b9649642a +Subproject commit 4fe67293f98407fb98a7a1cb14bc5ef55ed48443 diff --git a/scm/data/processed_case_input/gabls3.nc b/scm/data/processed_case_input/gabls3.nc index 70247b335..415e99fe8 100644 Binary files a/scm/data/processed_case_input/gabls3.nc and b/scm/data/processed_case_input/gabls3.nc differ diff --git a/scm/data/processed_case_input/gabls3_noahmp.nc b/scm/data/processed_case_input/gabls3_noahmp.nc index 8c09c72f1..615d2cfb5 100644 Binary files a/scm/data/processed_case_input/gabls3_noahmp.nc and b/scm/data/processed_case_input/gabls3_noahmp.nc differ diff --git a/scm/data/processed_case_input/gabls3_ruc.nc b/scm/data/processed_case_input/gabls3_ruc.nc new file mode 100644 index 000000000..4ca083fca Binary files /dev/null and b/scm/data/processed_case_input/gabls3_ruc.nc differ diff --git a/scm/etc/case_config/gabls3_ruc.nml b/scm/etc/case_config/gabls3_ruc.nml new file mode 100644 index 000000000..b3c0c1835 --- /dev/null +++ b/scm/etc/case_config/gabls3_ruc.nml @@ -0,0 +1,29 @@ +$case_config +model_name = 'FV3', +n_columns = 1, +case_name = 'gabls3_ruc', +dt = 600.0, +time_scheme = 1, +runtime = 86400, +n_itt_out = 1, +n_itt_diag = 6, +n_levels = 64, +output_file = 'output', +case_data_dir = '../data/processed_case_input', +vert_coord_data_dir = '../data/vert_coord_data', +thermo_forcing_type = 2, +mom_forcing_type = 2, +relax_time = 7200.0, +sfc_flux_spec = .false., +lsm_ics = .true., +do_spinup = .true., +spinup_timesteps = 12, +sfc_roughness_length_cm = 15.0, +sfc_type = 1, +reference_profile_choice = 2, +year = 2006, +month = 7, +day = 1, +hour = 12, +column_area = 1.45E8, +$end diff --git a/scm/etc/scm_qsub_example.py b/scm/etc/scm_qsub_example.py index e3b4edaca..457013dff 100755 --- a/scm/etc/scm_qsub_example.py +++ b/scm/etc/scm_qsub_example.py @@ -21,7 +21,7 @@ USER = os.getenv('USER') STRINGS = [USER, 'ucar.edu'] MY_EMAIL = '@'.join(STRINGS) -PROC = Popen('qsub', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=True) +PROC = Popen('qsub -V', shell=True, stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=True) ### Begin User-editable section ### JOB_NAME = "test_job" diff --git a/scm/etc/scripts/GABLS3_LSM.ncl b/scm/etc/scripts/GABLS3_LSM.ncl index ac16fc514..85f4dc2a9 100644 --- a/scm/etc/scripts/GABLS3_LSM.ncl +++ b/scm/etc/scripts/GABLS3_LSM.ncl @@ -431,7 +431,11 @@ begin tg3@units = "K" tg3@long_name = "deep soil temperature" - zorlw = 15.0 ;(cm) from case specs + zorl = 15.0 ;(cm) from case specs + zorl@long_name = "composite surface roughness length" + zorl@units = "cm" + + zorlw = zorl ;(cm) from case specs zorlw@long_name = "surface roughness length over ocean" zorlw@units = "cm" @@ -482,7 +486,7 @@ begin ;derive friction velocity from known values wind1 = sqrt(u_specified(1)^2 + v_specified(1)^2) - uustar = von_K*wind1/log(u_h(1)/(0.01*zorlw)) + uustar = von_K*wind1/log(u_h(1)/(0.01*zorl)) uustar@units = "m s-1" uustar@long_name = "friction velocity" @@ -533,11 +537,11 @@ begin tsfcl@long_name = "surface skin temperature over land" tsfcl@units = "K" - zorll = zorlw ;(cm) from case specs + zorll = zorl ;(cm) from case specs zorll@long_name = "surface roughness length over land" zorll@units = "cm" - zorli = zorlw ;(cm) from case specs + zorli = zorl ;(cm) from case specs zorli@long_name = "surface roughness length over ice" zorli@units = "cm" @@ -545,7 +549,7 @@ begin stc@units = "K" stc@long_name = "initial profile of soil temperature" - smc = (/ 0.197, 0.197, 0.197, 0.197 /) + smc = (/ 0.203, 0.203, 0.203, 0.203 /) smc@units = "m3 m-3" smc@long_name = "initial profile of soil moisture" @@ -636,6 +640,10 @@ begin filevarattdef(g1,"tg3",tg3) g1->tg3 = tg3 + filevardef(g1,"zorl",typeof(zorl),"ncl_scalar") + filevarattdef(g1,"zorl",zorl) + g1->zorl = zorl + filevardef(g1,"zorlw",typeof(zorlw),"ncl_scalar") filevarattdef(g1,"zorlw",zorlw) g1->zorlw = zorlw diff --git a/scm/etc/scripts/GABLS3_LSM_NoahMP.ncl b/scm/etc/scripts/GABLS3_LSM_NoahMP.ncl index 621c39bf1..aa956aebf 100644 --- a/scm/etc/scripts/GABLS3_LSM_NoahMP.ncl +++ b/scm/etc/scripts/GABLS3_LSM_NoahMP.ncl @@ -431,7 +431,11 @@ begin tg3@units = "K" tg3@long_name = "deep soil temperature" - zorlw = 15.0 ;(cm) from case specs + zorl = 15.0 ;(cm) from case specs + zorl@long_name = "composite surface roughness length" + zorl@units = "cm" + + zorlw = zorl ;(cm) from case specs zorlw@long_name = "surface roughness length over ocean" zorlw@units = "cm" @@ -482,7 +486,7 @@ begin ;derive friction velocity from known values wind1 = sqrt(u_specified(1)^2 + v_specified(1)^2) - uustar = von_K*wind1/log(u_h(1)/(0.01*zorlw)) + uustar = von_K*wind1/log(u_h(1)/(0.01*zorl)) uustar@units = "m s-1" uustar@long_name = "friction velocity" @@ -533,11 +537,11 @@ begin tsfcl@long_name = "surface skin temperature over land" tsfcl@units = "K" - zorll = zorlw ;(cm) from case specs + zorll = zorl ;(cm) from case specs zorll@long_name = "surface roughness length over land" zorll@units = "cm" - zorli = zorlw ;(cm) from case specs + zorli = zorl ;(cm) from case specs zorli@long_name = "surface roughness length over ice" zorli@units = "cm" @@ -636,6 +640,10 @@ begin filevarattdef(g1,"tg3",tg3) g1->tg3 = tg3 + filevardef(g1,"zorl",typeof(zorl),"ncl_scalar") + filevarattdef(g1,"zorl",zorl) + g1->zorl = zorl + filevardef(g1,"zorlw",typeof(zorlw),"ncl_scalar") filevarattdef(g1,"zorlw",zorlw) g1->zorlw = zorlw diff --git a/scm/etc/scripts/GABLS3_LSM_RUC.ncl b/scm/etc/scripts/GABLS3_LSM_RUC.ncl new file mode 100644 index 000000000..167cb5023 --- /dev/null +++ b/scm/etc/scripts/GABLS3_LSM_RUC.ncl @@ -0,0 +1,1018 @@ +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" + +begin + +;Define constants + + missing_value = -9999.0 + g = 9.80665 ;gravity (m s-2) + R_dry = 287. ;ideal gas dry air constant (J kg-1 K-1) + R_vap = 461.5 ;gas constant for water vapor + c_p = 1004. ;specific heat at const pressure (J kg-1 K-1) + l_v = 2.5e6 ;latent heat of vaporization (J kg-1) + l_s = 2.836e6 ;latent heat of sublimation (J kg-1) + P0 = 100000. ;intial pressure (Pa) + von_K = 0.4 ;vonKarmann constant + kappa = R_dry/c_p + +; GABLS initial vertical profiles and forcings. For levels with no entries, +; values are linearly interpolated based on nearby points +; Data source: http://projects.knmi.nl/gabls/setup.html + +; Time, 01Jul2006 12:00:00 UTC - 02Jul2006 12:00:00 UTC +; Must be in seconds + time = fspan(0,86400,145) + time@long_name = "elapsed time since the beginning of the simulation" + time@units = "s" + +; timb = fspan(0,86400,145) +; print(timb) + +; Lat and Lon of the column + lat = 51.9711 + lat@long_name = "latitude of column" + lat@units = "degrees N" + lon = 4.9267 + lon@long_name = "longitude of column" + lon@units = "degrees E" + +; Height levels for initial values of temperature and specific humidity +; Other variables will need to be interpolated + height = (/ 0., 2., 10., 20., 40., 80., 140., 205., 1800., 2200., \ + 5000., 12000., 14000., 50000. /) + height@long_name = "physical height at pressure levels" + height@units = "m" + +; Soil levels for RUC LSM + soil_depth = (/ 0.0, 0.01, 0.04, 0.10, 0.30, 0.60, 1.0, 1.6, 3.0 /) + soil_depth@long_name = "depth of bottom of soil layers" + soil_depth@units = "m" + +; Total water specific humidity (q) at heights listed above + qt = (/ 9.3e-3, 9.3e-3, 8.5e-3, 8.4e-3, 8.3e-3, 8.2e-3, 8.1e-3, \ + 8.0e-3, 7.5e-3, 2.0e-3, 0.3e-3, 0.01e-3, 0.003e-3, 0. /) + qt@long_name = "initial profile of total water specific humidity" + qt@units = "kg kg^-1" + +; Converting levels (m) to pressures (Pa) using Exner eqn +; Calculate Exner and pressure at surface + p_surf_hpa = 1024.4 ;(hPa, given on GABLS setup link) + p_surf = new((/ dimsizes(time) /),float) + p_surf(0) = p_surf_hpa * 100. ;now in Pa, rather than hPa + + do i=1,dimsizes(time)-1,1 + p_surf(i) = p_surf(0) + end do + + p_surf@long_name = "surface pressure" + p_surf@units = "Pa" + + levels = new((/ dimsizes(height) /),float) + levels@long_name = "pressure levels" + levels@units = "Pa" + +; Air temperature profile (not needed in netCDF input file, but necessary +; for calculating other variables that gmtb-scm will find necessary) + T_C = (/ 27.0, 27.0, 26.4, 26.2, 25.9, 25.5, 24.8, 24.3, 9.0, 9.0, \ + -6.4, -61., -54., -50. /) + T_K = T_C + 273.15 ;convert from degC to Kelvin + + T_surf = new((/ dimsizes(time) /),float) + T_surf = 23.4 + 273.15 + +; do i=0,dimsizes(time)-1,1 +; T_surf(i) = T_K(0) +; end do + T_surf@long_name = "surface absolute temperature" + T_surf@units = "K" + +; Calculate approximate virtual temperature to solve hypsometric equation +; for the model required pressure levels + T_V = new ((/ dimsizes(levels) /),float) + do i=0,dimsizes(levels)-1,1 + T_V(i) = T_K(i)*((1 + (qt(i)/0.622))/(1+qt(i))) + end do + +; Solve the hypsometric equation for p2 to attain each new pressure level +; and fill the array (the variable to fill is levels, created above.) +; The value for levels(0) = p_surf(0) + levels(0) = p_surf(0) + do i=1,dimsizes(levels)-1,1 + levels(i) = levels(i-1) / (exp((g*(height(i)-height(i-1)))/(R_dry \ + * T_V(i)))) + end do + +; Calculate potential temperature (theta) in order to calculate the +; necessary thetail (ice-liquid water potential temperature) + theta = new((/ dimsizes(levels) /),float) + do i=0,dimsizes(levels)-1,1 + theta(i) = T_K(i)*((P0/levels(i))^(kappa)) + end do + +; Calculate thetail. Because ql and qi are not given, this is not possible, +; so the initial profile of thetail will be assumed theta. +; The remainder of the equation used in the gmtb-scm user's guide would be +; zero, leaving you with only theta. + thetail = theta + thetail@long_name = "initial profile of ice-liquid water potential "+ \ + "temperature" + thetail@units = "K" + +; Horizontal motion (u and v winds), needs to be interpolated + u_specified = (/ 0.0, -4.0, -5.5, -5.5, -2.0, -2.0 /) + u_h = (/ 0., 10., 353., 1238., 2000., 50000. /) + + u = linint1(u_h,u_specified,False,height,0) + u@long_name = "initial profile of E-W horizontal wind" + u@units = "m s^-1" + + v_specified = (/ 0.0, -0.4, -0.5, -0.5, 2., 2. /) + v_h = u_h + + v = linint1(v_h,v_specified,False,height,0) + v@long_name = "initial profile of N-S horizontal wind" + v@units = "m s^-1" + +; Use GLDAS to derive the sensible and latent heat flux for the Cabauw +; location, then you can compare the GLDAS H/Le to the given intial Bowen +; ratio. Begin by reading in the data from each file. +; gl_dir = "/glade/scratch/damico/GLDAS_GABLS3/" +; gl_00 = addfile(gl_dir+"GLDAS_NOAH025_3H.A20060701.0000.021.nc4.SUB.nc4", \ +; "r") + +; gl_lat = gl_00->lat +; gl_lon = gl_00->lon + +; glny = dimsizes(gl_lat) +; glnx = dimsizes(gl_lon) + +; delete(gl_00) + +; start_date = 0. +; end_date = 0. + +; start_date@units = "hours since 2006-07-01 12:00:00" +; end_date@units = "hours since 2006-07-02 12:00:00" +; gldt = 3 + +; end_date_on_start = cd_convert(end_date,start_date@units) + +; gl_time_hr = ispan(0,toint(end_date_on_start),gldt) +; gl_time_hr@units = start_date@units + +; gl_time = tofloat(gl_time_hr) * 60 * 60 +; gl_time@units = "seconds since 2006-07-01 12:00:00" + +; time_yyyymmddhh = toint(cd_calendar(gl_time,0)) + +; glnt = dimsizes(gl_time_hr) + +; gl_shf = new((/ glnt,glny,glnx /),float) +; gl_lhf = new((/ glnt,glny,glnx /),float) +; do i=0,glnt-1,1 +; fname = gl_dir + "GLDAS_NOAH025_3H.A" + \ +; sprinti("%0.4i",time_yyyymmddhh(i,0)) + \ +; sprinti("%0.2i",time_yyyymmddhh(i,1)) + \ +; sprinti("%0.2i",time_yyyymmddhh(i,2)) +"."+ \ +; sprinti("%0.2i",time_yyyymmddhh(i,3))+ \ +; sprinti("%0.2i",time_yyyymmddhh(i,4))+ \ +; ".021.nc4.SUB.nc4" +; if (isfilepresent(fname)) then +; gl_ff = addfile(fname,"r") +; gl_shf(i,:,:) = gl_ff->Qh_tavg(0,:,:) +; gl_lhf(i,:,:) = gl_ff->Qle_tavg(0,:,:) +; print("Read->"+fname) +; else +; print(fname+" MISSING") +; end if +; end do + +; Interpolate GLDAS shf and lhf to the single Cabauw point. +; shf_pt = linint2_points_Wrap(gl_lon,gl_lat,gl_shf,False,lon,lat,0) +; lhf_pt = linint2_points_Wrap(gl_lon,gl_lat,gl_lhf,False,lon,lat,0) + +; Interpolate with time +; shf = linint1(gl_time,shf_pt(:,0),False,time,0) +; lhf = linint1(gl_time,lhf_pt(:,0),False,time,0) + +; Convert from W m^-2 to K m s^-1 and kg kg^-1 m s^-1 +; sh_flux_sfc = shf * R_dry * T_surf / (c_p * p_surf) +; lh_flux_sfc = lhf * R_dry * T_surf / (l_v * p_surf) + +; sh_flux_sfc@long_name = "surface sensible heat flux" +; sh_flux_sfc@units = "K m s^-1" +; lh_flux_sfc@long_name = "surface latent heat flux" +; lh_flux_sfc@units = "kg kg^-1 m s^-1" + +; Calculate vertical motion, needed in both w and omega, given in omega. Omega +; begins at 0.12 Pa s^-1 at 1500 m at the beginning of the study period. After +; five hours, omega starts to decrease to 0, ending at 0 after seven hours. +; Omega is measured at 1500 m, and decreases from that value to 0 m. It also +; decreases from 0.12 to 5000 m, above which the value is zero. + om_specified = (/ 0.12, 0.12, 0., 0. /) + omega_t = (/ 0,18000,25200, 86400 /) + om = linint1(omega_t,om_specified,False,time,0) + + omega = new((/ dimsizes(levels),dimsizes(time) /),float) + omega_h = (/ 0.,1500.,5000.,50000. /) + omega_int = new((/ dimsizes(omega_h),dimsizes(time) /),float) + do i=0,dimsizes(time)-1,1 + omega_int(:,i) = (/ 0.,om(i),0.,0. /) + omega(:,i) = linint1(omega_h,omega_int(:,i),False,height,0) + end do + + omega@long_name = "large scale pressure vertical velocity" + omega@units = "Pa s^-1" + +; Convert omega into w, also required input +; Begin by solving for rho + rho = new((/ dimsizes(levels) /),float) + do i=0,dimsizes(levels)-1,1 + rho(i) = levels(i)/(R_dry*T_K(i)) + end do + +; Solve for w (equation commented below) +; w = -om / (rho*g) + w_ls = new((/ dimsizes(levels),dimsizes(time) /),float) + do i=0,dimsizes(levels)-1,1 + do j=0,dimsizes(time)-1,1 + w_ls(i,j) = omega(i,j) / (rho(i)*g) + if (.not. ismissing(w_ls(i,j)) .and. w_ls(i,j) .ne. 0.) then + w_ls(i,j) = w_ls(i,j) * -1. + end if + end do + end do + + w_ls@long_name = "large scale vertical velocity" + w_ls@units = "m s^-1" + +; Geostrophic winds are given at certain times but only at the surface. +; Instructions from GABLS suggest interpolation to surface winds up to +; 2000 m, and then constant from that point. This should be done for each +; time step. Starting with time. + u_geo_spec = (/ -7.8, -7.8, -6.5, -5.0, -5.0, -6.5 /) + u_geo_t = (/ 0, 21600, 39600, 54000, 64800, 86400 /) + u_geo = linint1(u_geo_t,u_geo_spec,False,time,0) + + v_geo_spec = (/ 0.0, 0.0, 4.5, 4.5, 4.5, 2.5 /) + v_geo_t = (/ 0, 21600, 39600, 54000, 64800, 86400 /) + v_geo = linint1(v_geo_t,v_geo_spec,False,time,0) + +; Interpolate with height, linearly, from the surface values above to -2.0 +; m s^-1 for u and 2.0 m s^-1 for v. Above 2000 m, geostrophic wind is +; constant with height (according to GABLS guide). + u_g = new((/ dimsizes(levels),dimsizes(time) /),float) + u_geo_h = (/ 0., 2000., 50000. /) + u_geo_int = new((/ dimsizes(u_geo_h),dimsizes(time) /),float) + + do i=0,dimsizes(time)-1,1 + u_geo_int(:,i) = (/ u_geo(i), -2., -2. /) + u_g(:,i) = linint1(u_geo_h,u_geo_int(:,i),False,height,0) + end do + + v_g = new((/ dimsizes(levels),dimsizes(time) /),float) + v_geo_int = new((/ dimsizes(u_geo_h),dimsizes(time) /),float) + + v_geo_h = u_geo_h + + do i=0,dimsizes(time)-1,1 + v_geo_int(:,i) = (/ v_geo(i), 2., 2. /) + v_g(:,i) = linint1(v_geo_h,v_geo_int(:,i),False,height,0) + end do + + u_g@long_name = "large scale geostrophic E-W wind" + v_g@long_name = "large scale geostrophic N-S wind" + u_g@units = "m s^-1" + v_g@unit = "m s^-1" + +; Horizontal temperature tendency due to advection, which can be converted to +; thetail tendency. The values are give; between 200-1000m, decreasing +; to zero down towards the surface and decreasing to zero from 1000m to +; 1500m, with zero entirely above. Once again, there needs to be +; interpolation in time and with height. + T_advec_spec = (/ -2.5e-5,7.5e-5, 0., 0. /) + T_advec_time = (/ 0, 46800,64800,86400 /) + T_advec = linint1(T_advec_time,T_advec_spec,False,time,0) + +; Interpolate with height + h_advec_T = new((/ dimsizes(levels),dimsizes(time) /),float) + h_advec_h = (/ 0., 200., 1000., 1500., 50000. /) + T_advec_int = new((/ dimsizes(h_advec_h),dimsizes(time) /),float) + + do i=0,dimsizes(time)-1,1 + T_advec_int(:,i) = (/ 0.,T_advec(i),T_advec(i),0.,0. /) + h_advec_T(:,i) = linint1(h_advec_h,T_advec_int(:,i),False,height,0) + end do + +; Covert to theta_il + h_advec_thetail = new((/ dimsizes(levels),dimsizes(time) /),float) + do i=0,dimsizes(levels)-1,1 + h_advec_thetail(i,:) = h_advec_T(i,:) * ((P0 / levels(i))^(kappa)) + end do + + h_advec_thetail@long_name = "prescribed theta_il tendency due to "+ \ + "horizontal advection" + h_advec_thetail@units = "K s^-1" + +; The last remaining variable given in the data that the gmtb-scm can +; use is horizontal specific humidity tendency. The values are given +; between 200-1000m, decreasing to zero down towards the surface and +; decreasing to zero from 1000m to 1500, with zero entirely above. Once +; again, there needs to be interpolation in time and with height. + h_advec_spec = (/ 0., 8.e-8, 0., -8.e-8, 0., 0. /) + h_advec_time = (/ 0, 32400, 43200, 50400, 61200, 86400 /) + h_advec = linint1(h_advec_time,h_advec_spec,False,time,0) + +; Interpolate with height + h_advec_qt = new((/ dimsizes(levels),dimsizes(time) /),float) +; h_advec_h = (/ 0., 200., 1000., 1500., 50000. /) + h_advec_int = new((/ dimsizes(h_advec_h),dimsizes(time) /),float) + + do i=0,dimsizes(time)-1,1 + h_advec_int(:,i) = (/ 0., h_advec(i), h_advec(i), 0., 0. /) + h_advec_qt(:,i) = linint1(h_advec_h,h_advec_int(:,i),False,height,0) + end do + + h_advec_qt@long_name = "prescribed q_t tendency due to horizontal " + \ + "advection" + h_advec_qt@units = "kg kg^-1 s^-1" + +; The following variables are not in the GABLS data available or cannot +; be derived with the information given. They are still required by the +; model. +; Initial profiles + ql = new((/ dimsizes(levels) /),float) + ql = 0. + ql@long_name = "initial profile of liquid water specific humidity" + ql@units = "kg kg^-1" + + qi = new((/ dimsizes(levels) /),float) + qi = 0. + qi@long_name = "initial profile of ice water specific humidity" + qi@units = "kg kg^-1" + + tke = new((/ dimsizes(levels) /),float) + tke = 0. + tke@long_name = "initial profile of turbulence kinetic energy" + tke@units = "m^2 s^-2" + + ozone = new((/ dimsizes(levels) /),float) + ozone = 0. + ozone@long_name = "initial profile of ozone mass mixing ratio" + ozone@units = "kg kg^-1" + +; Forcing + + u_nudge = new((/ dimsizes(levels),dimsizes(time) /),float) + u_nudge = 0. + u_nudge@long_name = "E-W wind to nudge toward" + u_nudge@units = "m s^-1" + + v_nudge = new((/ dimsizes(levels),dimsizes(time) /),float) + v_nudge = 0. + v_nudge@long_name = "N-S wind to nudge toward" + v_nudge@units = "m s^-1" + + T_nudge = new((/ dimsizes(levels),dimsizes(time) /),float) + T_nudge = 0. + T_nudge@long_name = "absolute temperature to nudge toward" + T_nudge@units = "K" + + thil_nudge = new((/ dimsizes(levels),dimsizes(time) /),float) + thil_nudge = 0. + thil_nudge@long_name = "potential temperature to nudge toward" + thil_nudge@units = "K" + + qt_nudge = new((/ dimsizes(levels),dimsizes(time) /),float) + qt_nudge = 0. + qt_nudge@long_name = "q_t to nudge toward" + qt_nudge@units = "kg kg^-1" + + dT_dt_rad = new((/ dimsizes(levels),dimsizes(time) /),float) + dT_dt_rad = 0. + dT_dt_rad@long_name = "prescribed radiative heating rate" + dT_dt_rad@units = "K s^-1" + +; h_advec_thetail = new((/ dimsizes(levels),dimsizes(time) /), \ +; float) +; h_advec_thetail = 0. +; h_advec_thetail@long_name = "prescribed theta_il tendency due to " + \ +; "horizontal advection" +; h_advec_thetail@units = "K s^-1" + + v_advec_thetail = new((/ dimsizes(levels),dimsizes(time) /), \ + float) + v_advec_thetail = 0. + v_advec_thetail@long_name = "prescribed theta_il tendency due to " + \ + "vertical advection" + v_advec_thetail@units = "K s^-1" + + v_advec_qt = new((/ dimsizes(levels),dimsizes(time) /),float) + v_advec_qt = 0. + v_advec_qt@long_name = "prescribe q_t tendency due to vertical "+\ + "advection" + v_advec_qt@units = "kg kg^-1 s^-1" + +; Noah LSM initialization + slmsk = 1.0 ;corresponds to land from the UFS ICs + slmsk@long_name = "land-sea-ice mask" + + tsfco = 305.0 + tsfco@long_name = "sea surface temperature OR surface skin temperature over land OR sea ice surface skin temperature (depends on value of slmsk)" + tsfco@units = "K" + + weasd = 0.0 + weasd@long_name = "water equivalent accumulated snow depth" + weasd@units = "mm" + + tg3 = 283.15 ;corresponds to case specs + tg3@units = "K" + tg3@long_name = "deep soil temperature" + + zorl = 15.0 ;(cm) from case specs + zorl@long_name = "composite surface roughness length" + zorl@units = "cm" + + zorlw = zorl ;(cm) from case specs + zorlw@long_name = "surface roughness length over ocean" + zorlw@units = "cm" + + alvsf = 0.23 ;corresponds to case specs + alvsf@long_name = "60 degree vis albedo with strong cosz dependency" + + alnsf = 0.23 ;corresponds to case specs + alnsf@long_name = "60 degree nir albedo with strong cosz dependency" + + alvwf = 0.23 ;corresponds to case specs + alvwf@long_name = "60 degree vis albedo with weak cosz dependency" + + alnwf = 0.23 ;corresponds to case specs + alnwf@long_name = "60 degree nir albedo with weak cosz dependency" + + facsf = 0.50556319952 ;corresponds to value from UFS ICs (static from fix files for C768) + facsf@long_name = "fractional coverage with strong cosz dependency" + + facwf = 0.49443680048 ;corresponds to value from UFS ICs (static from fix files for C768) + facwf@long_name = "fractional coverage with weak cosz dependency" + + vegfrac = 0.75 ; the value of 100% specified in the case + ; specifications clashes with the maximum value + ; for this gridpoint from UFS ICs + vegfrac@long_name = "vegetation fraction" + + canopy = 0.5 ;corresponds to typical value for grassland in the growing season + canopy@units = "kg m-2" + canopy@long_name = "amount of water stored in canopy" + + f10m = missing_value ;this is not required + f10m@long_name = "ratio of sigma level 1 wind and 10m wind" + + t2m = missing_value; this is not required + t2m@long_name = "2-meter absolute temperature" + t2m@units = "K" + + q2m = missing_value; this is not required + q2m@long_name = "2-meter specific humidity" + q2m@units = "kg kg-1" + + vegtyp = 10 ;corresponds to the "grasslands" type for the IGBP (vegsrc=1) + ;dataset + vegtyp@long_name = "vegetation type (1-12)" + + soiltyp = 12 ;corresponds to "clay" soil type + soiltyp@long_name = "soil type (1-12)" + + ;derive friction velocity from known values + wind1 = sqrt(u_specified(1)^2 + v_specified(1)^2) + uustar = von_K*wind1/log(u_h(1)/(0.01*zorl)) + uustar@units = "m s-1" + uustar@long_name = "friction velocity" + + ffmm = missing_value ;this is not required + ffmm@long_name = "Monin-Obukhov similarity function for momentum" + + ffhh = missing_value ;this is not required + ffhh@long_name = "Monin-Obukhov similarity function for heat" + + hice = 0.0 ;corresponds to value from UFS ICs + hice@units = "m" + hice@long_name = "sea ice thickness" + + fice = 0.0 ;corresponds to value from UFS ICs + fice@long_name = "ice fraction" + + tisfc = tsfco ; there is no ice, but this corresponds to the surface skin temperature + tisfc@units = "K" + tisfc@long_name = "ice surface temperature" + + tprcp = missing_value; not required + tprcp@units = "m" + tprcp@long_name = "instantaneous total precipitation amount" + + srflag = missing_value; not required + srflag@long_name = "snow/rain flag for precipitation" + + snwdph = 0.0 ;corresponds to case specs + snwdph@units = "mm" + snwdph@long_name = "water equivalent snow depth" + + shdmin = 0.01 ; corresponds to minimum vegetation fraction from surface fixed data (is not actually used in Noah LSM) + shdmin@long_name = "minimum vegetation fraction" + + shdmax = 0.8 ; the maximum value from the surface fixed data is closer to 0.75 (is not actually used in Noah LSM) + shdmax@long_name = "maximum vegetation fraction" + + slopetyp = 1 + slopetyp@long_name = "slope type (1-9)" + + snoalb = 0.728796124458 ;corresponds to value from UFS ICs (static from fix files for C768) + snoalb@long_name = "maximum snow albedo" + + sncovr = 0.0 ;corresponds to case specs + sncovr@long_name = "snow area fraction" + + tsfcl = tsfco + tsfcl@long_name = "surface skin temperature over land" + tsfcl@units = "K" + + zorll = zorl ;(cm) from case specs + zorll@long_name = "surface roughness length over land" + zorll@units = "cm" + + zorli = zorl ;(cm) from case specs + zorli@long_name = "surface roughness length over ice" + zorli@units = "cm" + + wetness = 0.37; from rucinit using Noah ICs + wetness@long_name = "normalized soil wetness for land surface model" + + clw_surf_land = 0.0 + clw_surf_land@long_name = "moist cloud water mixing ratio at surface over land" + clw_surf_land@units = "kg kg-1" + + clw_surf_ice = 0.0 + clw_surf_ice@long_name = "moist cloud water mixing ratio at surface over ice" + clw_surf_ice@units = "kg kg-1" + + qwv_surf_land = qt(0)/(1.0 - qt(0)) + qwv_surf_land@long_name = "water vapor mixing ratio at surface over land" + qwv_surf_land@units = "kg kg-1" + + qwv_surf_ice = qwv_surf_land + qwv_surf_ice@long_name = "water vapor mixing ratio at surface over ice" + qwv_surf_ice@units = "kg kg-1" + + tsnow_land = 273.14 + tsnow_land@long_name = "snow temperature at the bottom of the first snow layer over land" + tsnow_land@units = "K" + + tsnow_ice = tsnow_land + tsnow_ice@long_name = "snow temperature at the bottom of the first snow layer over ice" + tsnow_ice@units = "K" + + snowfallac_land = 0.0 + snowfallac_land@long_name = "run-total snow accumulation on the ground" + snowfallac_land@units = "kg m-2" + + snowfallac_ice = 0.0 + snowfallac_ice@long_name = "run-total snow accumulation on the ice" + snowfallac_ice@units = "kg m-2" + + sncovr_ice = 0.0 + sncovr_ice@long_name = "surface snow area fraction over ice" + + sfalb_lnd = 0.23 + sfalb_lnd@long_name = "mean surface diffused sw albedo over land" + + sfalb_lnd_bck = 0.23 + sfalb_lnd_bck@long_name = "surface snow-free albedo over ice" + + sfalb_ice = 0.6 + sfalb_ice@long_name = "mean surface diffused sw albedo over ice" + + emis_ice = 0.97 + emis_ice@long_name = "surface lw emissivity in fraction over ice" + + tslb = (/ 296.55, 295.95, 294.55, 292.55, 290.65, 288.39, 285.35, 284.03, 283.15/) + tslb@units = "K" + tslb@long_name = "initial profile of soil temperature for RUC LSM" + + smois = (/ 0.31, 0.31, 0.31, 0.31, 0.31, 0.31, 0.31, 0.31, 0.31/) + smois@units = "fraction" + smois@long_name = "volumetric fraction of soil moisture for RUC LSM" + + sh2o = smois + sh2o@units = "fraction" + sh2o@long_name = "volume fraction of unfrozen soil moisture for RUC LSM" + + smfr = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) + smfr@units = "fraction" + smfr@long_name = "volume fraction of frozen soil moisture for RUC LSM" + + flfr = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) + flfr@long_name = "flag for frozen soil physics for RUC LSM" + +; Define dimension sizes for file creation + ntim = dimsizes(time) + klev = dimsizes(levels) + ksoil = dimsizes(soil_depth) + +; Write the netCDF output file, titled "gabls3_ruc.nc" + setfileoption("nc","FileStructure","Advanced") + setfileoption("nc","Format","NetCDF4") + fout = "gabls3_ruc.nc" ;File name of output file + system("/bin/rm -fv "+fout) ;If a version already exists, delete it. + fo = addfile(fout,"c") ;Create the netCDF file + + setfileoption(fo,"DefineMode",True) ;Enter define mode + fAtt = True ;Set file attributes to True + fAtt@title = "GMTB SCM forcing file for GABLS3 case" + fAtt@creation_date = systemfunc("date") + fileattdef(fo,fAtt) ;Add attibutes to output file + + dimNames = (/ "time","levels","nsoil" /) + dimSizes = (/ ntim,klev,ksoil /) + dimUnlim = (/ False,True,True /) +; dimNames = (/ "levels","time" /) +; dimSizes = (/ klev,ntim /) +; dimUnlim = (/ True,False /) + + filedimdef(fo,dimNames,dimSizes,dimUnlim) + +; chunkSizes = (/ dimsizes(time),dimsizes(levels)/2 /) +; filechunkdimdef(fo,dimNames,chunkSizes,dimUnlim) + + filevardef(fo,"time",typeof(time),"time") + filevarattdef(fo,"time",time) + fo->time = (/ time /) + + filevardef(fo,"levels",typeof(levels),"levels") + filevarattdef(fo,"levels",levels) + fo->levels = (/ levels /) + + filevardef(fo,"soil_depth",typeof(soil_depth),"nsoil") + filevarattdef(fo,"soil_depth",soil_depth) + fo->soil_depth = (/ soil_depth /) + +; Define group names for netCDF4 file + grpnames = (/ "scalars","initial","forcing" /) + filegrpdef(fo,grpnames) + +; Scalars --> Lat and Lon + + g1 = fo=>/scalars + g2 = fo=>/initial + g3 = fo=>/forcing + +; filedimdef(g1,dimNames,dimSizes,dimUnlim) + +; Initial + + filedimdef(g2,dimNames,dimSizes,dimUnlim) + filedimdef(g3,dimNames,dimSizes,dimUnlim) + + filevardef(g1,"lat",typeof(lat),"ncl_scalar") + filevarattdef(g1,"lat",lat) + g1->lat = lat + + filevardef(g1,"lon",typeof(lon),"ncl_scalar") + filevarattdef(g1,"lon",lon) + g1->lon = lon + + filevardef(g1,"slmsk",typeof(slmsk),"ncl_scalar") + filevarattdef(g1,"slmsk",slmsk) + g1->slmsk = slmsk + + filevardef(g1,"tsfco",typeof(tsfco),"ncl_scalar") + filevarattdef(g1,"tsfco",tsfco) + g1->tsfco = tsfco + + filevardef(g1,"weasd",typeof(weasd),"ncl_scalar") + filevarattdef(g1,"weasd",weasd) + g1->weasd = weasd + + filevardef(g1,"tg3",typeof(tg3),"ncl_scalar") + filevarattdef(g1,"tg3",tg3) + g1->tg3 = tg3 + + filevardef(g1,"zorl",typeof(zorl),"ncl_scalar") + filevarattdef(g1,"zorl",zorl) + g1->zorl = zorl + + filevardef(g1,"zorlw",typeof(zorlw),"ncl_scalar") + filevarattdef(g1,"zorlw",zorlw) + g1->zorlw = zorlw + + filevardef(g1,"alvsf",typeof(alvsf),"ncl_scalar") + filevarattdef(g1,"alvsf",alvsf) + g1->alvsf = alvsf + + filevardef(g1,"alnsf",typeof(alnsf),"ncl_scalar") + filevarattdef(g1,"alnsf",alnsf) + g1->alnsf = alnsf + + filevardef(g1,"alvwf",typeof(alvwf),"ncl_scalar") + filevarattdef(g1,"alvwf",alvwf) + g1->alvwf = alvwf + + filevardef(g1,"alnwf",typeof(alnwf),"ncl_scalar") + filevarattdef(g1,"alnwf",alnwf) + g1->alnwf = alnwf + + filevardef(g1,"facsf",typeof(facsf),"ncl_scalar") + filevarattdef(g1,"facsf",facsf) + g1->facsf = facsf + + filevardef(g1,"facwf",typeof(facwf),"ncl_scalar") + filevarattdef(g1,"facwf",facwf) + g1->facwf = facwf + + filevardef(g1,"vegfrac",typeof(vegfrac),"ncl_scalar") + filevarattdef(g1,"vegfrac",vegfrac) + g1->vegfrac = vegfrac + + filevardef(g1,"canopy",typeof(canopy),"ncl_scalar") + filevarattdef(g1,"canopy",canopy) + g1->canopy = canopy + + filevardef(g1,"f10m",typeof(f10m),"ncl_scalar") + filevarattdef(g1,"f10m",f10m) + g1->f10m = f10m + + filevardef(g1,"t2m",typeof(t2m),"ncl_scalar") + filevarattdef(g1,"t2m",t2m) + g1->t2m = t2m + + filevardef(g1,"q2m",typeof(q2m),"ncl_scalar") + filevarattdef(g1,"q2m",q2m) + g1->q2m = q2m + + filevardef(g1,"vegtyp",typeof(vegtyp),"ncl_scalar") + filevarattdef(g1,"vegtyp",vegtyp) + g1->vegtyp = vegtyp + + filevardef(g1,"soiltyp",typeof(soiltyp),"ncl_scalar") + filevarattdef(g1,"soiltyp",soiltyp) + g1->soiltyp = soiltyp + + filevardef(g1,"uustar",typeof(uustar),"ncl_scalar") + filevarattdef(g1,"uustar",uustar) + g1->uustar = uustar + + filevardef(g1,"ffmm",typeof(ffmm),"ncl_scalar") + filevarattdef(g1,"ffmm",ffmm) + g1->ffmm = ffmm + + filevardef(g1,"ffhh",typeof(ffhh),"ncl_scalar") + filevarattdef(g1,"ffhh",ffhh) + g1->ffhh = ffhh + + filevardef(g1,"hice",typeof(hice),"ncl_scalar") + filevarattdef(g1,"hice",hice) + g1->hice = hice + + filevardef(g1,"fice",typeof(fice),"ncl_scalar") + filevarattdef(g1,"fice",fice) + g1->fice = fice + + filevardef(g1,"tisfc",typeof(tisfc),"ncl_scalar") + filevarattdef(g1,"tisfc",tisfc) + g1->tisfc = tisfc + + filevardef(g1,"tprcp",typeof(tprcp),"ncl_scalar") + filevarattdef(g1,"tprcp",tprcp) + g1->tprcp = tprcp + + filevardef(g1,"srflag",typeof(srflag),"ncl_scalar") + filevarattdef(g1,"srflag",srflag) + g1->srflag = srflag + + filevardef(g1,"snwdph",typeof(snwdph),"ncl_scalar") + filevarattdef(g1,"snwdph",snwdph) + g1->snwdph = snwdph + + filevardef(g1,"shdmin",typeof(shdmin),"ncl_scalar") + filevarattdef(g1,"shdmin",shdmin) + g1->shdmin = shdmin + + filevardef(g1,"shdmax",typeof(shdmax),"ncl_scalar") + filevarattdef(g1,"shdmax",shdmax) + g1->shdmax = shdmax + + filevardef(g1,"slopetyp",typeof(slopetyp),"ncl_scalar") + filevarattdef(g1,"slopetyp",slopetyp) + g1->slopetyp = slopetyp + + filevardef(g1,"snoalb",typeof(snoalb),"ncl_scalar") + filevarattdef(g1,"snoalb",snoalb) + g1->snoalb = snoalb + + filevardef(g1,"sncovr",typeof(sncovr),"ncl_scalar") + filevarattdef(g1,"sncovr",sncovr) + g1->sncovr = sncovr + + filevardef(g1,"tsfcl",typeof(tsfcl),"ncl_scalar") + filevarattdef(g1,"tsfcl",tsfcl) + g1->tsfcl = tsfcl + + filevardef(g1,"zorll",typeof(zorll),"ncl_scalar") + filevarattdef(g1,"zorll",zorll) + g1->zorll = zorll + + filevardef(g1,"zorli",typeof(zorli),"ncl_scalar") + filevarattdef(g1,"zorli",zorli) + g1->zorli = zorli + + filevardef(g1,"wetness",typeof(wetness),"ncl_scalar") + filevarattdef(g1,"wetness",wetness) + g1->wetness = wetness + + filevardef(g1,"clw_surf_land",typeof(clw_surf_land),"ncl_scalar") + filevarattdef(g1,"clw_surf_land",clw_surf_land) + g1->clw_surf_land = clw_surf_land + + filevardef(g1,"clw_surf_ice",typeof(clw_surf_ice),"ncl_scalar") + filevarattdef(g1,"clw_surf_ice",clw_surf_ice) + g1->clw_surf_ice = clw_surf_ice + + filevardef(g1,"qwv_surf_land",typeof(qwv_surf_land),"ncl_scalar") + filevarattdef(g1,"qwv_surf_land",qwv_surf_land) + g1->qwv_surf_land = qwv_surf_land + + filevardef(g1,"qwv_surf_ice",typeof(qwv_surf_ice),"ncl_scalar") + filevarattdef(g1,"qwv_surf_ice",qwv_surf_ice) + g1->qwv_surf_ice = qwv_surf_ice + + filevardef(g1,"tsnow_land",typeof(tsnow_land),"ncl_scalar") + filevarattdef(g1,"tsnow_land",tsnow_land) + g1->tsnow_land = tsnow_land + + filevardef(g1,"tsnow_ice",typeof(tsnow_ice),"ncl_scalar") + filevarattdef(g1,"tsnow_ice",tsnow_ice) + g1->tsnow_ice = tsnow_ice + + filevardef(g1,"snowfallac_land",typeof(snowfallac_land),"ncl_scalar") + filevarattdef(g1,"snowfallac_land",snowfallac_land) + g1->snowfallac_land = snowfallac_land + + filevardef(g1,"snowfallac_ice",typeof(snowfallac_ice),"ncl_scalar") + filevarattdef(g1,"snowfallac_ice",snowfallac_ice) + g1->snowfallac_ice = snowfallac_ice + + filevardef(g1,"sncovr_ice",typeof(sncovr_ice),"ncl_scalar") + filevarattdef(g1,"sncovr_ice",sncovr_ice) + g1->sncovr_ice = sncovr_ice + + filevardef(g1,"sfalb_lnd",typeof(sfalb_lnd),"ncl_scalar") + filevarattdef(g1,"sfalb_lnd",sfalb_lnd) + g1->sfalb_lnd = sfalb_lnd + + filevardef(g1,"sfalb_lnd_bck",typeof(sfalb_lnd_bck),"ncl_scalar") + filevarattdef(g1,"sfalb_lnd_bck",sfalb_lnd_bck) + g1->sfalb_lnd_bck = sfalb_lnd_bck + + filevardef(g1,"sfalb_ice",typeof(sfalb_ice),"ncl_scalar") + filevarattdef(g1,"sfalb_ice",sfalb_ice) + g1->sfalb_ice = sfalb_ice + + filevardef(g1,"emis_ice",typeof(emis_ice),"ncl_scalar") + filevarattdef(g1,"emis_ice",emis_ice) + g1->emis_ice = emis_ice + + filevardef(g2,"height",typeof(height),"levels") + filevarattdef(g2,"height",height) + g2->height = (/height/) + + filevardef(g2,"thetail",typeof(thetail),"levels") + filevarattdef(g2,"thetail",thetail) + g2->thetail = (/thetail/) + + filevardef(g2,"qt",typeof(qt),"levels") + filevarattdef(g2,"qt",qt) + g2->qt = (/qt/) + + filevardef(g2,"ql",typeof(ql),"levels") + filevarattdef(g2,"ql",ql) + g2->ql = (/ql/) + + filevardef(g2,"qi",typeof(qi),"levels") + filevarattdef(g2,"qi",qi) + g2->qi = (/qi/) + + filevardef(g2,"u",typeof(u),"levels") + filevarattdef(g2,"u",u) + g2->u = (/u/) + + filevardef(g2,"v",typeof(v),"levels") + filevarattdef(g2,"v",v) + g2->v = (/v/) + + filevardef(g2,"tke",typeof(tke),"levels") + filevarattdef(g2,"tke",tke) + g2->tke = (/tke/) + + filevardef(g2,"ozone",typeof(ozone),"levels") + filevarattdef(g2,"ozone",ozone) + g2->ozone = (/ozone/) + + filevardef(g2,"tslb",typeof(tslb),"nsoil") + filevarattdef(g2,"tslb",tslb) + g2->tslb = (/tslb/) + + filevardef(g2,"smois",typeof(smois),"nsoil") + filevarattdef(g2,"smois",smois) + g2->smois = (/smois/) + + filevardef(g2,"sh2o",typeof(sh2o),"nsoil") + filevarattdef(g2,"sh2o",sh2o) + g2->sh2o = (/sh2o/) + + filevardef(g2,"smfr",typeof(smfr),"nsoil") + filevarattdef(g2,"smfr",smfr) + g2->smfr = (/smfr/) + + filevardef(g2,"flfr",typeof(flfr),"nsoil") + filevarattdef(g2,"flfr",flfr) + g2->flfr = (/flfr/) + +; Forcing + + filevardef(g3,"p_surf",typeof(p_surf),"time") + filevarattdef(g3,"p_surf",p_surf) + g3->p_surf = (/ p_surf /) + + filevardef(g3,"T_surf",typeof(T_surf),"time") + filevarattdef(g3,"T_surf",T_surf) + g3->T_surf = (/ T_surf /) + +; filevardef(g3,"sh_flux_sfc",typeof(sh_flux_sfc),"time") +; filevarattdef(g3,"sh_flux_sfc",sh_flux_sfc) +; g3->sh_flux_sfc = (/ sh_flux_sfc /) + +; filevardef(g3,"lh_flux_sfc",typeof(lh_flux_sfc),"time") +; filevarattdef(g3,"lh_flux_sfc",lh_flux_sfc) +; g3->lh_flux_sfc = (/ lh_flux_sfc /) + + filevardef(g3,"w_ls",typeof(w_ls),(/"levels","time"/)) + filevarattdef(g3,"w_ls",w_ls) + g3->w_ls = (/ w_ls /) + + filevardef(g3,"omega",typeof(omega),(/"levels","time"/)) + filevarattdef(g3,"omega",omega) + g3->omega = (/ omega /) + + filevardef(g3,"u_g",typeof(u_g),(/"levels","time"/)) + filevarattdef(g3,"u_g",u_g) + g3->u_g = (/ u_g /) + + filevardef(g3,"v_g",typeof(v_g),(/"levels","time"/)) + filevarattdef(g3,"v_g",v_g) + g3->v_g = (/ v_g /) + + filevardef(g3,"u_nudge",typeof(u_nudge),(/"levels","time"/)) + filevarattdef(g3,"u_nudge",u_nudge) + g3->u_nudge = (/ u_nudge /) + + filevardef(g3,"v_nudge",typeof(v_nudge),(/"levels","time"/)) + filevarattdef(g3,"v_nudge",v_nudge) + g3->v_nudge = (/ v_nudge /) + + filevardef(g3,"T_nudge",typeof(T_nudge),(/"levels","time"/)) + filevarattdef(g3,"T_nudge",T_nudge) + g3->T_nudge = (/ T_nudge /) + + filevardef(g3,"thil_nudge",typeof(thil_nudge),(/"levels","time"/)) + filevarattdef(g3,"thil_nudge",thil_nudge) + g3->thil_nudge = (/ thil_nudge /) + + filevardef(g3,"qt_nudge",typeof(qt_nudge),(/"levels","time"/)) + filevarattdef(g3,"qt_nudge",qt_nudge) + g3->qt_nudge = (/ qt_nudge /) + + filevardef(g3,"dT_dt_rad",typeof(dT_dt_rad),(/"levels","time"/)) + filevarattdef(g3,"dT_dt_rad",dT_dt_rad) + g3->dT_dt_rad = (/ dT_dt_rad /) + + filevardef(g3,"h_advec_thetail",typeof(h_advec_thetail), \ + (/"levels","time"/)) + filevarattdef(g3,"h_advec_thetail",h_advec_thetail) + g3->h_advec_thetail = (/ h_advec_thetail /) + + filevardef(g3,"v_advec_thetail",typeof(v_advec_thetail), \ + (/"levels","time"/)) + filevarattdef(g3,"v_advec_thetail",v_advec_thetail) + g3->v_advec_thetail = (/ v_advec_thetail /) + + filevardef(g3,"h_advec_qt",typeof(h_advec_qt),(/"levels","time"/)) + filevarattdef(g3,"h_advec_qt",h_advec_qt) + g3->h_advec_qt = (/ h_advec_qt /) + + filevardef(g3,"v_advec_qt",typeof(v_advec_qt),(/"levels","time"/)) + filevarattdef(g3,"v_advec_qt",v_advec_qt) + g3->v_advec_qt = (/ v_advec_qt /) + + print("wrote new file: "+fout) + +end diff --git a/scm/src/GFS_typedefs.F90 b/scm/src/GFS_typedefs.F90 index c2acb3c80..b50183da9 100644 --- a/scm/src/GFS_typedefs.F90 +++ b/scm/src/GFS_typedefs.F90 @@ -1287,6 +1287,7 @@ module GFS_typedefs integer :: kdt !< current forecast iteration logical :: first_time_step !< flag signaling first time step for time integration routine logical :: restart !< flag whether this is a coldstart (.false.) or a warmstart/restart (.true.) + logical :: lsm_cold_start logical :: hydrostatic !< flag whether this is a hydrostatic or non-hydrostatic run integer :: jdat(1:8) !< current forecast date and time !< (yr, mon, day, t-zone, hr, min, sec, mil-sec) @@ -4811,6 +4812,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%kdt = nint(Model%fhour*con_hr/Model%dtp) Model%first_time_step = .true. Model%restart = restart + Model%lsm_cold_start = .not. restart Model%hydrostatic = hydrostatic Model%jdat(1:8) = jdat(1:8) allocate(Model%si(Model%levs+1)) @@ -5956,6 +5958,7 @@ subroutine control_print(Model) print *, ' sec : ', Model%sec print *, ' first_time_step : ', Model%first_time_step print *, ' restart : ', Model%restart + print *, ' lsm_cold_start : ', Model%lsm_cold_start print *, ' hydrostatic : ', Model%hydrostatic endif diff --git a/scm/src/GFS_typedefs.meta b/scm/src/GFS_typedefs.meta index 63ec41b69..7900294f8 100644 --- a/scm/src/GFS_typedefs.meta +++ b/scm/src/GFS_typedefs.meta @@ -5237,6 +5237,12 @@ units = flag dimensions = () type = logical +[lsm_cold_start] + standard_name = do_lsm_cold_start + long_name = flag to signify LSM is cold-started + units = flag + dimensions = () + type = logical [hydrostatic] standard_name = flag_for_hydrostatic_solver long_name = flag for hydrostatic solver from dynamics diff --git a/scm/src/scm_input.F90 b/scm/src/scm_input.F90 index b05e0ac54..648ba0243 100644 --- a/scm/src/scm_input.F90 +++ b/scm/src/scm_input.F90 @@ -219,86 +219,73 @@ subroutine get_case_init(scm_state, scm_input) type(scm_state_type), intent(in) :: scm_state type(scm_input_type), target, intent(inout) :: scm_input - integer :: input_nlev !< number of levels in the input file - integer :: input_nsoil !< number of soil levels in the input file - integer :: input_ntimes !< number of times represented in the input file - integer :: input_nsnow !< number of snow levels in the input file - integer :: input_nice !< number of sea ice levels in the input file - integer :: input_nsoil_plus_nsnow !< number of combined snow and soil levels in the input file - - ! dimension variables - real(kind=dp), allocatable :: input_pres(:) !< input file pressure levels (Pa) - real(kind=dp), allocatable :: input_time(:) !< input file times (seconds since the beginning of the case) - - !initial profile variables - real(kind=dp), allocatable :: input_thetail(:) !< ice-liquid water potential temperature profile (K) - real(kind=dp), allocatable :: input_temp(:) !< temperature profile (K) - real(kind=dp), allocatable :: input_qt(:) !< total water specific humidity profile (kg kg^-1) - real(kind=dp), allocatable :: input_ql(:) !< liquid water specific humidity profile (kg kg^-1) - real(kind=dp), allocatable :: input_qi(:) !< ice water specific humidity profile (kg kg^-1) - real(kind=dp), allocatable :: input_u(:) !< east-west horizontal wind profile (m s^-1) - real(kind=dp), allocatable :: input_v(:) !< north-south horizontal wind profile (m s^-1) - real(kind=dp), allocatable :: input_tke(:) !< TKE profile (m^2 s^-2) - real(kind=dp), allocatable :: input_ozone(:) !< ozone profile (kg kg^-1) -! additional land info - real(kind=dp), allocatable :: input_stc(:) !< soil temperature (K) - real(kind=dp), allocatable :: input_smc(:) !< total soil moisture content (fraction) - real(kind=dp), allocatable :: input_slc(:) !< liquid soil moisture content (fraction) - !real(kind=dp), allocatable :: input_pres_i(:) !< interface pressures - !real(kind=dp), allocatable :: input_pres_l(:) !< layer pressures - real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) - real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) - real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) - real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) - real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) - real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) - real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) - real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) - real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) - real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) - real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics - - integer :: input_vegsrc !< vegetation source - integer :: input_vegtyp !< vegetation type - integer :: input_soiltyp!< soil type - integer :: input_slopetype !< slope type + integer :: input_nlev !< number of levels in the input file + integer :: input_nsoil !< number of soil levels in the input file + integer :: input_nsnow !< number of snow levels in the input file + integer :: input_nice !< number of sea ice levels in the input file + integer :: input_ntimes !< number of times represented in the input file + integer :: input_nsoil_plus_nsnow !< number of combined snow and soil levels in the input file real(kind=dp) :: input_lat !< column latitude (deg) real(kind=dp) :: input_lon !< column longitude (deg) - real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) - real(kind=dp) :: input_vegfrac !< vegetation fraction - real(kind=dp) :: input_shdmin !< minimun vegetation fraction - real(kind=dp) :: input_shdmax !< maximun vegetation fraction - real(kind=dp) :: input_zorlw !< surfce roughness length over water [cm] - real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] - real(kind=dp) :: input_canopy !< amount of water stored in canopy (kg m-2) - real(kind=dp) :: input_hice !< sea ice thickness (m) - real(kind=dp) :: input_fice !< ice fraction (frac) - real(kind=dp) :: input_tisfc !< ice surface temperature (K) - real(kind=dp) :: input_snwdph !< water equivalent snow depth (mm) - real(kind=dp) :: input_snoalb !< maximum snow albedo (frac) - real(kind=dp) :: input_sncovr !< snow area fraction (frac) real(kind=dp) :: input_area !< surface area [m^2] + + integer :: input_vegsrc !< vegetation source + + real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] + real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) + real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) real(kind=dp) :: input_tg3 !< deep soil temperature (K) - real(kind=dp) :: input_uustar !< surface friction velocity (m s-1) + real(kind=dp) :: input_zorl !< composite surface roughness length (cm) real(kind=dp) :: input_alvsf !< 60 degree vis albedo with strong cosz dependency - real(kind=dp) :: input_alnsf !< 60 degree nir albedo with strong cosz dependency real(kind=dp) :: input_alvwf !< 60 degree vis albedo with weak cosz dependency + real(kind=dp) :: input_alnsf !< 60 degree nir albedo with strong cosz dependency real(kind=dp) :: input_alnwf !< 60 degree nir albedo with weak cosz dependency real(kind=dp) :: input_facsf !< fractional coverage with strong cosz dependency real(kind=dp) :: input_facwf !< fractional coverage with weak cosz dependency - real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) + real(kind=dp) :: input_vegfrac !< vegetation fraction + real(kind=dp) :: input_canopy !< amount of water stored in canopy (kg m-2) real(kind=dp) :: input_f10m !< ratio of sigma level 1 wind and 10m wind real(kind=dp) :: input_t2m !< 2-meter absolute temperature (K) real(kind=dp) :: input_q2m !< 2-meter specific humidity (kg kg-1) + integer :: input_vegtyp !< vegetation type + integer :: input_soiltyp!< soil type + real(kind=dp) :: input_uustar !< surface friction velocity (m s-1) real(kind=dp) :: input_ffmm !< Monin-Obukhov similarity function for momentum real(kind=dp) :: input_ffhh !< Monin-Obukhov similarity function for heat + real(kind=dp) :: input_hice !< sea ice thickness (m) + real(kind=dp) :: input_fice !< ice fraction (frac) + real(kind=dp) :: input_tisfc !< ice surface temperature (K) real(kind=dp) :: input_tprcp !< instantaneous total precipitation amount (m) real(kind=dp) :: input_srflag !< snow/rain flag for precipitation + real(kind=dp) :: input_snwdph !< water equivalent snow depth (mm) + real(kind=dp) :: input_shdmin !< minimun vegetation fraction + real(kind=dp) :: input_shdmax !< maximun vegetation fraction + integer :: input_slopetype !< slope type + real(kind=dp) :: input_snoalb !< maximum snow albedo (frac) + real(kind=dp) :: input_sncovr !< snow area fraction (frac) + real(kind=dp) :: input_snodl !< snowd on land portion of cell + real(kind=dp) :: input_weasdl !< weasd on land portion of cell + real(kind=dp) :: input_tsfc !< tsfc composite real(kind=dp) :: input_tsfcl !< surface skin temperature over land (K) + real(kind=dp) :: input_zorlw !< surfce roughness length over water [cm] real(kind=dp) :: input_zorll !< surface roughness length over land (cm) real(kind=dp) :: input_zorli !< surface roughness length over ice (cm) + real(kind=dp) :: input_albdirvis_lnd !< + real(kind=dp) :: input_albdirnir_lnd !< + real(kind=dp) :: input_albdifvis_lnd !< + real(kind=dp) :: input_albdifnir_lnd !< + real(kind=dp) :: input_emis_lnd !< + real(kind=dp) :: input_albdirvis_ice !< + real(kind=dp) :: input_albdirnir_ice !< + real(kind=dp) :: input_albdifvis_ice !< + real(kind=dp) :: input_albdifnir_ice !< real(kind=dp) :: input_zorlwav !< surface roughness length from wave model (cm) + real(kind=dp), allocatable :: input_stc(:) !< soil temperature (K) + real(kind=dp), allocatable :: input_smc(:) !< total soil moisture content (fraction) + real(kind=dp), allocatable :: input_slc(:) !< liquid soil moisture content (fraction) + real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) + real(kind=dp) :: input_stddev !< standard deviation of subgrid orography (m) real(kind=dp) :: input_convexity !< convexity of subgrid orography real(kind=dp) :: input_ol1 !< fraction of grid box with subgrid orography higher than critical height 1 @@ -318,7 +305,7 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_landfrac !< fraction of horizontal grid area occupied by land real(kind=dp) :: input_lakefrac !< fraction of horizontal grid area occupied by lake real(kind=dp) :: input_lakedepth !< lake depth (m) - + real(kind=dp) :: input_tvxy !< vegetation temperature (K) real(kind=dp) :: input_tgxy !< ground temperature for Noahmp (K) real(kind=dp) :: input_tahxy !< canopy air temperature (K) @@ -348,11 +335,12 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_deeprechxy !< recharge to or from the water table when deep (m) real(kind=dp) :: input_rechxy !< recharge to or from the water table when shallow (m) real(kind=dp) :: input_snowxy !< number of snow layers - real(kind=dp) :: input_albdvis !< - real(kind=dp) :: input_albdnir !< - real(kind=dp) :: input_albivis !< - real(kind=dp) :: input_albinir !< - real(kind=dp) :: input_emiss !< + + real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) + real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) + real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) + real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) + real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) real(kind=dp) :: input_tref !< sea surface reference temperature for NSST (K) real(kind=dp) :: input_z_c !< sub-layer cooling thickness for NSST (m) @@ -383,15 +371,38 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_snowfallac_land !< run-total snow accumulation on the ground over land for RUC LSM (kg m-2) real(kind=dp) :: input_snowfallac_ice !< run-total snow accumulation on the ground over ice for RUC LSM (kg m-2) real(kind=dp) :: input_sncovr_ice !< + real(kind=dp) :: input_sfalb_lnd + real(kind=dp) :: input_sfalb_lnd_bck + real(kind=dp) :: input_sfalb_ice + real(kind=dp) :: input_emis_ice real(kind=dp) :: input_lai !< leaf area index for RUC LSM - !surface time-series variables + real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) + real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics + + ! dimension variables + !real(kind=dp), allocatable :: input_pres_i(:) !< interface pressures + !real(kind=dp), allocatable :: input_pres_l(:) !< layer pressures + real(kind=dp), allocatable :: input_pres(:) !< input file pressure levels (Pa) + real(kind=dp), allocatable :: input_time(:) !< input file times (seconds since the beginning of the case) + + !initial profile variables + real(kind=dp), allocatable :: input_thetail(:) !< ice-liquid water potential temperature profile (K) + real(kind=dp), allocatable :: input_temp(:) !< temperature profile (K) + real(kind=dp), allocatable :: input_qt(:) !< total water specific humidity profile (kg kg^-1) + real(kind=dp), allocatable :: input_ql(:) !< liquid water specific humidity profile (kg kg^-1) + real(kind=dp), allocatable :: input_qi(:) !< ice water specific humidity profile (kg kg^-1) + real(kind=dp), allocatable :: input_u(:) !< east-west horizontal wind profile (m s^-1) + real(kind=dp), allocatable :: input_v(:) !< north-south horizontal wind profile (m s^-1) + real(kind=dp), allocatable :: input_tke(:) !< TKE profile (m^2 s^-2) + real(kind=dp), allocatable :: input_ozone(:) !< ozone profile (kg kg^-1) + real(kind=dp), allocatable :: input_pres_surf(:) !< time-series of surface pressure (Pa) real(kind=dp), allocatable :: input_T_surf(:) !< time-series of surface temperature (K) - real(kind=dp), allocatable :: input_sh_flux_sfc(:) !< time-series of surface sensible heat flux (K m s^-1) - real(kind=dp), allocatable :: input_lh_flux_sfc(:) !< time-series of surface latent heat flux (kg kg^-1 m s^-1) - !2D (time, pressure) variables real(kind=dp), allocatable :: input_w_ls(:,:) !< 2D vertical velocity (m s^-1) real(kind=dp), allocatable :: input_omega(:,:) !< 2D pressure vertical velocity (Pa s^-1) real(kind=dp), allocatable :: input_u_g(:,:) !< 2D geostrophic east-west wind (m s^-1) @@ -407,6 +418,9 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp), allocatable :: input_v_advec_thetail(:,:) !< 2D theta_il tendency due to large-scale vertical advection (K s^-1) real(kind=dp), allocatable :: input_v_advec_qt(:,:) !< 2D q_t tendency due to large-scale horizontal vertical (kg kg^-1 s^-1) + real(kind=dp), allocatable :: input_sh_flux_sfc(:) !< time-series of surface sensible heat flux (K m s^-1) + real(kind=dp), allocatable :: input_lh_flux_sfc(:) !< time-series of surface latent heat flux (kg kg^-1 m s^-1) + CHARACTER(LEN=nf90_max_name) :: tmpName integer :: ncid, varID, grp_ncid, allocate_status,ierr real(kind=dp) :: nc_missing_value @@ -525,42 +539,55 @@ subroutine get_case_init(scm_state, scm_input) !possible scalars !Noah LSM parameters (when running with model ICs) - call NetCDF_read_var(grp_ncid, "vegsrc", .False., input_vegsrc ) - call NetCDF_read_var(grp_ncid, "vegtyp", .False., input_vegtyp ) - call NetCDF_read_var(grp_ncid, "soiltyp", .False., input_soiltyp ) - call NetCDF_read_var(grp_ncid, "slopetyp", .False., input_slopetype) - call NetCDF_read_var(grp_ncid, "tsfco", .False., input_tsfco) - call NetCDF_read_var(grp_ncid, "vegfrac", .False., input_vegfrac) - call NetCDF_read_var(grp_ncid, "shdmin", .False., input_shdmin) - call NetCDF_read_var(grp_ncid, "shdmax", .False., input_shdmax) - call NetCDF_read_var(grp_ncid, "zorlw", .False., input_zorlw) + call NetCDF_read_var(grp_ncid, "vegsrc", .False., input_vegsrc ) call NetCDF_read_var(grp_ncid, "slmsk", .False., input_slmsk) - call NetCDF_read_var(grp_ncid, "canopy", .False., input_canopy) - call NetCDF_read_var(grp_ncid, "hice", .False., input_hice) - call NetCDF_read_var(grp_ncid, "fice", .False., input_fice) - call NetCDF_read_var(grp_ncid, "tisfc", .False., input_tisfc) - call NetCDF_read_var(grp_ncid, "snwdph", .False., input_snwdph) - call NetCDF_read_var(grp_ncid, "snoalb", .False., input_snoalb) + call NetCDF_read_var(grp_ncid, "tsfco", .False., input_tsfco) + call NetCDF_read_var(grp_ncid, "weasd", .False., input_weasd) call NetCDF_read_var(grp_ncid, "tg3", .False., input_tg3) - call NetCDF_read_var(grp_ncid, "uustar", .False., input_uustar) + call NetCDF_read_var(grp_ncid, "zorl", .False., input_zorl) call NetCDF_read_var(grp_ncid, "alvsf", .False., input_alvsf) - call NetCDF_read_var(grp_ncid, "alnsf", .False., input_alnsf) call NetCDF_read_var(grp_ncid, "alvwf", .False., input_alvwf) + call NetCDF_read_var(grp_ncid, "alnsf", .False., input_alnsf) call NetCDF_read_var(grp_ncid, "alnwf", .False., input_alnwf) call NetCDF_read_var(grp_ncid, "facsf", .False., input_facsf) call NetCDF_read_var(grp_ncid, "facwf", .False., input_facwf) - call NetCDF_read_var(grp_ncid, "weasd", .False., input_weasd) + call NetCDF_read_var(grp_ncid, "vegfrac", .False., input_vegfrac) + call NetCDF_read_var(grp_ncid, "canopy", .False., input_canopy) call NetCDF_read_var(grp_ncid, "f10m", .False., input_f10m) call NetCDF_read_var(grp_ncid, "t2m", .False., input_t2m) call NetCDF_read_var(grp_ncid, "q2m", .False., input_q2m) + call NetCDF_read_var(grp_ncid, "vegtyp", .False., input_vegtyp ) + call NetCDF_read_var(grp_ncid, "soiltyp", .False., input_soiltyp ) + call NetCDF_read_var(grp_ncid, "uustar", .False., input_uustar) call NetCDF_read_var(grp_ncid, "ffmm", .False., input_ffmm) call NetCDF_read_var(grp_ncid, "ffhh", .False., input_ffhh) + call NetCDF_read_var(grp_ncid, "hice", .False., input_hice) + call NetCDF_read_var(grp_ncid, "fice", .False., input_fice) + call NetCDF_read_var(grp_ncid, "tisfc", .False., input_tisfc) call NetCDF_read_var(grp_ncid, "tprcp", .False., input_tprcp) call NetCDF_read_var(grp_ncid, "srflag", .False., input_srflag) + call NetCDF_read_var(grp_ncid, "snwdph", .False., input_snwdph) + call NetCDF_read_var(grp_ncid, "shdmin", .False., input_shdmin) + call NetCDF_read_var(grp_ncid, "shdmax", .False., input_shdmax) + call NetCDF_read_var(grp_ncid, "slopetyp",.False., input_slopetype) + call NetCDF_read_var(grp_ncid, "snoalb", .False., input_snoalb) call NetCDF_read_var(grp_ncid, "sncovr", .False., input_sncovr) + call NetCDF_read_var(grp_ncid, "snodl", .False., input_snodl) + call NetCDF_read_var(grp_ncid, "weasdl", .False., input_weasdl) + call NetCDF_read_var(grp_ncid, "tsfc", .False., input_tsfc) call NetCDF_read_var(grp_ncid, "tsfcl", .False., input_tsfcl) + call NetCDF_read_var(grp_ncid, "zorlw", .False., input_zorlw) call NetCDF_read_var(grp_ncid, "zorll", .False., input_zorll) call NetCDF_read_var(grp_ncid, "zorli", .False., input_zorli) + call NetCDF_read_var(grp_ncid, "albdirvis_lnd", .False., input_albdirvis_lnd) + call NetCDF_read_var(grp_ncid, "albdirnir_lnd", .False., input_albdirnir_lnd) + call NetCDF_read_var(grp_ncid, "albdifvis_lnd", .False., input_albdifvis_lnd) + call NetCDF_read_var(grp_ncid, "albdifnir_lnd", .False., input_albdifnir_lnd) + call NetCDF_read_var(grp_ncid, "emis_lnd", .False., input_emis_lnd) + call NetCDF_read_var(grp_ncid, "albdirvis_ice", .False., input_albdirvis_ice) + call NetCDF_read_var(grp_ncid, "albdirnir_ice", .False., input_albdirnir_ice) + call NetCDF_read_var(grp_ncid, "albdifvis_ice", .False., input_albdifvis_ice) + call NetCDF_read_var(grp_ncid, "albdifnir_ice", .False., input_albdifnir_ice) call NetCDF_read_var(grp_ncid, "zorlwav", .False., input_zorlwav) !orographic parameters @@ -614,11 +641,6 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "deeprechxy",.False., input_deeprechxy) call NetCDF_read_var(grp_ncid, "rechxy", .False., input_rechxy) call NetCDF_read_var(grp_ncid, "snowxy", .False., input_snowxy) - call NetCDF_read_var(grp_ncid, "albdvis", .False., input_albdvis) - call NetCDF_read_var(grp_ncid, "albdnir", .False., input_albdnir) - call NetCDF_read_var(grp_ncid, "albivis", .False., input_albivis) - call NetCDF_read_var(grp_ncid, "albinir", .False., input_albinir) - call NetCDF_read_var(grp_ncid, "emiss", .False., input_emiss) !NSST variables call NetCDF_read_var(grp_ncid, "tref", .False., input_tref) @@ -648,11 +670,14 @@ subroutine get_case_init(scm_state, scm_input) call NetCDF_read_var(grp_ncid, "qwv_surf_ice", .False., input_qwv_surf_ice) call NetCDF_read_var(grp_ncid, "tsnow_land", .False., input_tsnow_land) call NetCDF_read_var(grp_ncid, "tsnow_ice", .False., input_tsnow_ice) - call NetCDF_read_var(grp_ncid, "snowfall_acc_land",.False., input_snowfallac_land) - call NetCDF_read_var(grp_ncid, "snowfall_acc_ice", .False., input_snowfallac_ice) + call NetCDF_read_var(grp_ncid, "snowfallac_land", .False., input_snowfallac_land) + call NetCDF_read_var(grp_ncid, "snowfallac_ice", .False., input_snowfallac_ice) call NetCDF_read_var(grp_ncid, "sncovr_ice", .False., input_sncovr_ice) + call NetCDF_read_var(grp_ncid, "sfalb_lnd", .False., input_sfalb_lnd) + call NetCDF_read_var(grp_ncid, "sfalb_lnd_bck", .False., input_sfalb_lnd_bck) + call NetCDF_read_var(grp_ncid, "emis_ice", .False., input_emis_ice) call NetCDF_read_var(grp_ncid, "lai", .False., input_lai) - + !> - Read in the forcing data. !> - Find group ncid for forcing group. @@ -749,42 +774,55 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_flfr = input_flfr scm_input%input_vegsrc = input_vegsrc - scm_input%input_vegtyp = input_vegtyp - scm_input%input_soiltyp = input_soiltyp - scm_input%input_slopetype = input_slopetype - scm_input%input_tsfco = input_tsfco - scm_input%input_vegfrac = input_vegfrac - scm_input%input_shdmin = input_shdmin - scm_input%input_shdmax = input_shdmax - scm_input%input_zorlw = input_zorlw + scm_input%input_slmsk = input_slmsk - scm_input%input_canopy = input_canopy - scm_input%input_hice = input_hice - scm_input%input_fice = input_fice - scm_input%input_tisfc = input_tisfc - scm_input%input_snwdph = input_snwdph - scm_input%input_snoalb = input_snoalb - scm_input%input_sncovr = input_sncovr - scm_input%input_area = input_area + scm_input%input_tsfco = input_tsfco + scm_input%input_weasd = input_weasd scm_input%input_tg3 = input_tg3 - scm_input%input_uustar = input_uustar + scm_input%input_zorl = input_zorl scm_input%input_alvsf = input_alvsf - scm_input%input_alnsf = input_alnsf scm_input%input_alvwf = input_alvwf + scm_input%input_alnsf = input_alnsf scm_input%input_alnwf = input_alnwf scm_input%input_facsf = input_facsf scm_input%input_facwf = input_facwf - scm_input%input_weasd = input_weasd + scm_input%input_vegfrac = input_vegfrac + scm_input%input_canopy = input_canopy scm_input%input_f10m = input_f10m scm_input%input_t2m = input_t2m scm_input%input_q2m = input_q2m + scm_input%input_vegtyp = input_vegtyp + scm_input%input_soiltyp = input_soiltyp + scm_input%input_uustar = input_uustar scm_input%input_ffmm = input_ffmm scm_input%input_ffhh = input_ffhh + scm_input%input_hice = input_hice + scm_input%input_fice = input_fice + scm_input%input_tisfc = input_tisfc scm_input%input_tprcp = input_tprcp scm_input%input_srflag = input_srflag + scm_input%input_snwdph = input_snwdph + scm_input%input_shdmin = input_shdmin + scm_input%input_shdmax = input_shdmax + scm_input%input_slopetype = input_slopetype + scm_input%input_snoalb = input_snoalb + scm_input%input_sncovr = input_sncovr + scm_input%input_snodl = input_snodl + scm_input%input_weasdl = input_weasdl + scm_input%input_tsfc = input_tsfc scm_input%input_tsfcl = input_tsfcl + scm_input%input_zorlw = input_zorlw scm_input%input_zorll = input_zorll scm_input%input_zorli = input_zorli + scm_input%input_albdirvis_lnd = input_albdirvis_lnd + scm_input%input_albdirnir_lnd = input_albdirnir_lnd + scm_input%input_albdifvis_lnd = input_albdifvis_lnd + scm_input%input_albdifnir_lnd = input_albdifnir_lnd + scm_input%input_emis_lnd = input_emis_lnd + scm_input%input_albdirvis_ice = input_albdirvis_ice + scm_input%input_albdirnir_ice = input_albdirnir_ice + scm_input%input_albdifvis_ice = input_albdifvis_ice + scm_input%input_albdifnir_ice = input_albdifnir_ice scm_input%input_zorlwav = input_zorlwav scm_input%input_stddev = input_stddev @@ -836,11 +874,6 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_deeprechxy = input_deeprechxy scm_input%input_rechxy = input_rechxy scm_input%input_snowxy = input_snowxy - scm_input%input_albdvis = input_albdvis - scm_input%input_albdnir = input_albdnir - scm_input%input_albivis = input_albivis - scm_input%input_albinir = input_albinir - scm_input%input_emiss = input_emiss scm_input%input_tref = input_tref scm_input%input_z_c = input_z_c @@ -861,6 +894,8 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_dt_cool = input_dt_cool scm_input%input_qrain = input_qrain + scm_input%input_area = input_area + scm_input%input_wetness = input_wetness scm_input%input_clw_surf_land = input_clw_surf_land scm_input%input_clw_surf_ice = input_clw_surf_ice @@ -871,6 +906,9 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_snowfallac_land = input_snowfallac_land scm_input%input_snowfallac_ice = input_snowfallac_ice scm_input%input_sncovr_ice = input_sncovr_ice + scm_input%input_sfalb_lnd = input_sfalb_lnd + scm_input%input_sfalb_lnd_bck = input_sfalb_lnd_bck + scm_input%input_emis_ice = input_emis_ice scm_input%input_lai = input_lai !> @} diff --git a/scm/src/scm_output.F90 b/scm/src/scm_output.F90 index ce7e318d0..13c0feeec 100644 --- a/scm/src/scm_output.F90 +++ b/scm/src/scm_output.F90 @@ -91,7 +91,7 @@ subroutine output_init(scm_state, physics) CALL CHECK(NF90_DEF_DIM(NCID=ncid,NAME="vert_dim_layer",LEN=scm_state%n_levels,DIMID=vert_dim_id)) CALL CHECK(NF90_DEF_DIM(NCID=ncid,NAME="vert_dim_interface",LEN=scm_state%n_levels+1,DIMID=vert_dim_i_id)) CALL CHECK(NF90_DEF_DIM(NCID=ncid,NAME="vert_dim_rad",LEN=physics%Interstitial%lmk,DIMID=vert_dim_rad_id)) - CALL CHECK(NF90_DEF_DIM(NCID=ncid,NAME="vert_dim_soil",LEN=physics%Model%lsoil,DIMID=vert_dim_soil_id)) + CALL CHECK(NF90_DEF_DIM(NCID=ncid,NAME="vert_dim_soil",LEN=physics%Model%lsoil_lsm,DIMID=vert_dim_soil_id)) !> - Define the dimension variables. call NetCDF_def_var(ncid, 'time_inst', NF90_FLOAT, "model elapsed time for instantaneous variables", "s", dummy_id, (/ time_inst_id /)) @@ -225,10 +225,14 @@ subroutine output_init_sfcprop(ncid, time_inst_id, hor_dim_id, vert_dim_soil_id, integer :: i, dummy_id if (scm_state%model_ics .or. scm_state%lsm_ics) then - if (physics%Model%lsoil > 0) then - call NetCDF_def_var(ncid, 'soil_T', NF90_FLOAT, "soil temperature profile ", "K", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) - call NetCDF_def_var(ncid, 'soil_moisture', NF90_FLOAT, "soil moisture profile ", "m3 m-3", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) - call NetCDF_def_var(ncid, 'soil_moisture_unfrozen', NF90_FLOAT, "unfrozen soil moisture profile ", "m3 m-3", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) + if (physics%Model%lsoil_lsm > 0) then + call NetCDF_def_var(ncid, 'soil_T', NF90_FLOAT, "soil temperature profile ", "K", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) + call NetCDF_def_var(ncid, 'soil_moisture', NF90_FLOAT, "soil moisture profile ", "m3 m-3", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) + call NetCDF_def_var(ncid, 'soil_moisture_unfrozen', NF90_FLOAT, "unfrozen soil moisture profile ", "m3 m-3", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) + + call NetCDF_def_var(ncid, 'soil_moisture_total_vol_frac', NF90_FLOAT, "total soil moisture volume fraction profile ", "", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) + call NetCDF_def_var(ncid, 'soil_moisture_unfrozen_vol_frac', NF90_FLOAT, "unfrozen soil moisture volume fraction profile ", "", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) + call NetCDF_def_var(ncid, 'soil_moisture_frozen_vol_frac', NF90_FLOAT, "frozen soil moisture volume fraction profile ", "", dummy_id, (/ hor_dim_id, vert_dim_soil_id, time_inst_id /)) end if end if @@ -529,6 +533,7 @@ end subroutine output_append_forcing subroutine output_append_sfcprop(ncid, scm_state, physics) use scm_type_defs, only: scm_state_type, physics_type use NetCDF_put, only: NetCDF_put_var + use NetCDF_read, only: missing_value use scm_physical_constants, only: con_rd, con_fvirt, con_cp, con_hvap integer, intent(in) :: ncid @@ -537,12 +542,26 @@ subroutine output_append_sfcprop(ncid, scm_state, physics) integer :: i real(kind=dp), dimension(scm_state%n_cols) :: rho, shf, lhf + real(kind=dp), dimension(scm_state%n_cols,physics%Model%lsoil_lsm) :: missing_value_2D if (scm_state%model_ics .or. scm_state%lsm_ics) then - if (physics%Model%lsoil > 0) then - call NetCDF_put_var(ncid, 'soil_T', physics%Sfcprop%stc(:,:), scm_state%itt_out) - call NetCDF_put_var(ncid, 'soil_moisture', physics%Sfcprop%smc(:,:), scm_state%itt_out) - call NetCDF_put_var(ncid, 'soil_moisture_unfrozen', physics%Sfcprop%slc(:,:), scm_state%itt_out) + if (physics%Model%lsoil_lsm > 0) then + missing_value_2D = missing_value + if (physics%Model%lsm == physics%Model%lsm_noah .or. physics%Model%lsm == physics%Model%lsm_noahmp) then + call NetCDF_put_var(ncid, 'soil_T', physics%Sfcprop%stc(:,:), scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture', physics%Sfcprop%smc(:,:), scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_unfrozen', physics%Sfcprop%slc(:,:), scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_total_vol_frac', missing_value_2D, scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_unfrozen_vol_frac', missing_value_2D, scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_frozen_vol_frac', missing_value_2D, scm_state%itt_out) + else if (physics%Model%lsm == physics%Model%lsm_ruc) then + call NetCDF_put_var(ncid, 'soil_T', physics%Sfcprop%tslb(:,:), scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture', missing_value_2D, scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_unfrozen', missing_value_2D, scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_total_vol_frac', physics%Sfcprop%smois(:,:), scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_unfrozen_vol_frac', physics%Sfcprop%sh2o(:,:), scm_state%itt_out) + call NetCDF_put_var(ncid, 'soil_moisture_frozen_vol_frac', physics%Sfcprop%keepsmfr(:,:), scm_state%itt_out) + end if end if end if diff --git a/scm/src/scm_type_defs.F90 b/scm/src/scm_type_defs.F90 index 45f4beb86..dd333ac45 100644 --- a/scm/src/scm_type_defs.F90 +++ b/scm/src/scm_type_defs.F90 @@ -174,47 +174,67 @@ module scm_type_defs integer :: input_nsnow !< number of snow layers in the input file integer :: input_nice !< number of sea ice layers in the input file integer :: input_ntimes !< number of times in the input file where forcing is available - integer :: input_vegsrc !< vegetation source - integer :: input_vegtyp !< vegetation type classification - integer :: input_soiltyp !< - integer :: input_slopetype !< surface slope classification real(kind=dp) :: input_lat !< latitude of column center real(kind=dp) :: input_lon !< longitude of column center - real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) - real(kind=dp) :: input_vegfrac !< vegetation area fraction - real(kind=dp) :: input_shdmin !< minimun vegetation fraction - real(kind=dp) :: input_shdmax !< maximun vegetation fraction - real(kind=dp) :: input_zorlw !< surfce roughness length over water [cm] - real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] - real(kind=dp) :: input_canopy !< amount of water stored in canopy (kg m-2) - real(kind=dp) :: input_hice !< sea ice thickness (m) - real(kind=dp) :: input_fice !< ice fraction (frac) - real(kind=dp) :: input_tisfc !< ice surface temperature (K) - real(kind=dp) :: input_snwdph !< water equivalent snow depth (mm) - real(kind=dp) :: input_snoalb !< maximum snow albedo (frac) - real(kind=dp) :: input_sncovr !< snow area fraction (frac) real(kind=dp) :: input_area !< surface area [m^2] + + integer :: input_vegsrc !< vegetation source + + real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] + real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) + real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) real(kind=dp) :: input_tg3 !< deep soil temperature (K) - real(kind=dp) :: input_uustar !< surface friction velocity (m s-1) + real(kind=dp) :: input_zorl !< composite surface roughness length (cm) real(kind=dp) :: input_alvsf !< 60 degree vis albedo with strong cosz dependency - real(kind=dp) :: input_alnsf !< 60 degree nir albedo with strong cosz dependency real(kind=dp) :: input_alvwf !< 60 degree vis albedo with weak cosz dependency + real(kind=dp) :: input_alnsf !< 60 degree nir albedo with strong cosz dependency real(kind=dp) :: input_alnwf !< 60 degree nir albedo with weak cosz dependency real(kind=dp) :: input_facsf !< fractional coverage with strong cosz dependency real(kind=dp) :: input_facwf !< fractional coverage with weak cosz dependency - real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) + real(kind=dp) :: input_vegfrac !< vegetation area fraction + real(kind=dp) :: input_canopy !< amount of water stored in canopy (kg m-2) real(kind=dp) :: input_f10m !< ratio of sigma level 1 wind and 10m wind real(kind=dp) :: input_t2m !< 2-meter absolute temperature (K) real(kind=dp) :: input_q2m !< 2-meter specific humidity (kg kg-1) + integer :: input_vegtyp !< vegetation type classification + integer :: input_soiltyp !< + real(kind=dp) :: input_uustar !< surface friction velocity (m s-1) real(kind=dp) :: input_ffmm !< Monin-Obukhov similarity function for momentum real(kind=dp) :: input_ffhh !< Monin-Obukhov similarity function for heat + real(kind=dp) :: input_hice !< sea ice thickness (m) + real(kind=dp) :: input_fice !< ice fraction (frac) + real(kind=dp) :: input_tisfc !< ice surface temperature (K) real(kind=dp) :: input_tprcp !< instantaneous total precipitation amount (m) real(kind=dp) :: input_srflag !< snow/rain flag for precipitation + real(kind=dp) :: input_snwdph !< water equivalent snow depth (mm) + real(kind=dp) :: input_shdmin !< minimun vegetation fraction + real(kind=dp) :: input_shdmax !< maximun vegetation fraction + integer :: input_slopetype !< surface slope classification + real(kind=dp) :: input_snoalb !< maximum snow albedo (frac) + real(kind=dp) :: input_sncovr !< snow area fraction (frac) + real(kind=dp) :: input_snodl !< snowd on land portion of cell + real(kind=dp) :: input_weasdl !< weasd on land portion of cell + real(kind=dp) :: input_tsfc !< tsfc composite real(kind=dp) :: input_tsfcl !< surface skin temperature over land (K) + real(kind=dp) :: input_zorlw !< surfce roughness length over water [cm] real(kind=dp) :: input_zorll !< surface roughness length over land (cm) real(kind=dp) :: input_zorli !< surface roughness length over ice (cm) + real(kind=dp) :: input_albdirvis_lnd !< + real(kind=dp) :: input_albdirnir_lnd !< + real(kind=dp) :: input_albdifvis_lnd !< + real(kind=dp) :: input_albdifnir_lnd !< + real(kind=dp) :: input_emis_lnd !< + real(kind=dp) :: input_albdirvis_ice !< + real(kind=dp) :: input_albdifvis_ice !< + real(kind=dp) :: input_albdirnir_ice !< + real(kind=dp) :: input_albdifnir_ice !< real(kind=dp) :: input_zorlwav !< surface roughness length from wave model (cm) + real(kind=dp), allocatable :: input_stc(:) !< soil temperature (k) (initial) + real(kind=dp), allocatable :: input_smc(:) !< soil moisture content (g/g) (initial) + real(kind=dp), allocatable :: input_slc(:) !< soil liquid content (g/g) (initial) + real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) + real(kind=dp) :: input_stddev !< standard deviation of subgrid orography (m) real(kind=dp) :: input_convexity !< convexity of subgrid orography real(kind=dp) :: input_ol1 !< fraction of grid box with subgrid orography higher than critical height 1 @@ -264,11 +284,12 @@ module scm_type_defs real(kind=dp) :: input_deeprechxy !< recharge to or from the water table when deep (m) real(kind=dp) :: input_rechxy !< recharge to or from the water table when shallow (m) real(kind=dp) :: input_snowxy !< number of snow layers - real(kind=dp) :: input_albdvis !< - real(kind=dp) :: input_albdnir !< - real(kind=dp) :: input_albivis !< - real(kind=dp) :: input_albinir !< - real(kind=dp) :: input_emiss !< + + real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) + real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) + real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) + real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) + real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) real(kind=dp) :: input_tref !< sea surface reference temperature for NSST (K) real(kind=dp) :: input_z_c !< sub-layer cooling thickness for NSST (m) @@ -298,9 +319,19 @@ module scm_type_defs real(kind=dp) :: input_tsnow_ice !< snow temperature at the bottom of the first snow layer over ice for RUC LSM (K) real(kind=dp) :: input_snowfallac_land !< run-total snow accumulation on the ground over land for RUC LSM (kg m-2) real(kind=dp) :: input_snowfallac_ice !< run-total snow accumulation on the ground over ice for RUC LSM (kg m-2) - real(kind=dp) :: input_sncovr_ice !< + real(kind=dp) :: input_sncovr_ice !< + real(kind=dp) :: input_sfalb_lnd !< + real(kind=dp) :: input_sfalb_lnd_bck !< + real(kind=dp) :: input_sfalb_ice !< + real(kind=dp) :: input_emis_ice !< real(kind=dp) :: input_lai !< leaf area index for RUC LSM + real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) + real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics + real(kind=dp), allocatable :: input_pres_i(:) !< pressure (Pa) of input interface real(kind=dp), allocatable :: input_pres_l(:) !< pressure (Pa) of input levels real(kind=dp), allocatable :: input_pres(:) !< pressure (Pa) of input levels @@ -316,20 +347,6 @@ module scm_type_defs real(kind=dp), allocatable :: input_v(:) !< meridional wind (m/s) (initial) real(kind=dp), allocatable :: input_tke(:) !< turbulence kinetic energy (m^2/s^2) (initial) real(kind=dp), allocatable :: input_ozone(:) !< ozone mass mixing ratio (kg/kg) (initial) - real(kind=dp), allocatable :: input_stc(:) !< soil temperature (k) (initial) - real(kind=dp), allocatable :: input_smc(:) !< soil moisture conteng (g/g) (initial) - real(kind=dp), allocatable :: input_slc(:) !< soil liquid content (g/g) (initial) - real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) - real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) - real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) - real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) - real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) - real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) - real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) - real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) - real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) - real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) - real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics real(kind=dp), allocatable :: input_pres_surf(:) !< time-series of surface pressure (Pa) real(kind=dp), allocatable :: input_T_surf(:) !< time-series of surface temperture real(kind=dp), allocatable :: input_w_ls(:,:) !< large-scale vertical velocity (m/s) (time, levels) @@ -636,132 +653,90 @@ end subroutine scm_state_create subroutine scm_input_create(scm_input, ntimes, nlev, nsoil, nsnow, nice) class(scm_input_type) :: scm_input integer, intent(in) :: ntimes, nlev, nsoil, nsnow, nice - - scm_input%input_nlev = nlev + + scm_input%input_nlev = nlev + scm_input%input_nsoil = nsoil + scm_input%input_nsnow = nsnow + scm_input%input_nice = nice scm_input%input_ntimes = ntimes - scm_input%input_nsoil = nsoil - scm_input%input_nsnow = nsnow - scm_input%input_nice = nice - - allocate(scm_input%input_pres(nlev),scm_input%input_time(ntimes)) - scm_input%input_pres = real_zero - scm_input%input_time = real_zero - - allocate(scm_input%input_height(nlev), scm_input%input_thetail(nlev), scm_input%input_qt(nlev), scm_input%input_qv(nlev), & - scm_input%input_ql(nlev), scm_input%input_qi(nlev), scm_input%input_u(nlev), scm_input%input_v(nlev), scm_input%input_tke(nlev), & - scm_input%input_ozone(nlev), scm_input%input_temp(nlev)) - scm_input%input_height = real_zero - scm_input%input_thetail = real_zero - scm_input%input_temp = real_zero - scm_input%input_qt = real_zero - scm_input%input_qv = real_zero - scm_input%input_ql = real_zero - scm_input%input_qi = real_zero - scm_input%input_u = real_zero - scm_input%input_v = real_zero - scm_input%input_tke = real_zero - scm_input%input_ozone = real_zero - allocate(scm_input%input_pres_i(nlev+1),scm_input%input_pres_l(nlev)) - scm_input%input_pres_i = real_zero - scm_input%input_pres_l = real_zero - - allocate(scm_input%input_pres_surf(ntimes), & - scm_input%input_T_surf(ntimes), scm_input%input_sh_flux_sfc_kin(ntimes), scm_input%input_lh_flux_sfc_kin(ntimes), & - scm_input%input_sh_flux_sfc(ntimes), scm_input%input_lh_flux_sfc(ntimes), & - scm_input%input_w_ls(ntimes, nlev), scm_input%input_omega(ntimes, nlev), scm_input%input_u_g(ntimes, nlev), & - scm_input%input_v_g(ntimes, nlev), scm_input%input_dT_dt_rad(ntimes, nlev), scm_input%input_h_advec_thetail(ntimes, nlev), & - scm_input%input_h_advec_qt(ntimes, nlev), scm_input%input_v_advec_thetail(ntimes, nlev), & - scm_input%input_v_advec_qt(ntimes, nlev), scm_input%input_u_nudge(ntimes, nlev), scm_input%input_v_nudge(ntimes, nlev), & - scm_input%input_T_nudge(ntimes, nlev), scm_input%input_thil_nudge(ntimes, nlev), scm_input%input_qt_nudge(ntimes, nlev), & - scm_input%input_k_T_nudge(ntimes), scm_input%input_k_thil_nudge(ntimes), scm_input%input_k_qt_nudge(ntimes), & - scm_input%input_k_u_nudge(ntimes), scm_input%input_k_v_nudge(ntimes)) - scm_input%input_pres_surf = real_zero - scm_input%input_T_surf = real_zero - scm_input%input_lat = real_zero - scm_input%input_lon = real_zero + scm_input%input_lat = real_zero + scm_input%input_lon = real_zero + scm_input%input_area = real_zero - allocate(scm_input%input_tot_advec_t(ntimes, nlev), scm_input%input_tot_advec_theta(ntimes, nlev), & - scm_input%input_tot_advec_thetal(ntimes, nlev), scm_input%input_tot_advec_qv(ntimes, nlev), & - scm_input%input_tot_advec_u(ntimes, nlev), scm_input%input_tot_advec_v(ntimes, nlev), & - scm_input%input_pres_forcing(ntimes, nlev)) + scm_input%input_vegsrc = int_zero + + scm_input%input_slmsk = real_zero + scm_input%input_tsfco = real_zero + scm_input%input_weasd = real_zero + scm_input%input_tg3 = real_zero + scm_input%input_zorl = real_zero + scm_input%input_alvsf = real_zero + scm_input%input_alnsf = real_zero + scm_input%input_alvwf = real_zero + scm_input%input_alnwf = real_zero + scm_input%input_facsf = real_zero + scm_input%input_facwf = real_zero + scm_input%input_vegfrac = real_zero + scm_input%input_canopy = real_zero + scm_input%input_f10m = real_zero + scm_input%input_t2m = real_zero + scm_input%input_q2m = real_zero + scm_input%input_vegtyp = int_zero + scm_input%input_soiltyp = int_zero + scm_input%input_uustar = real_zero + scm_input%input_ffmm = real_zero + scm_input%input_ffhh = real_zero + scm_input%input_hice = real_zero + scm_input%input_fice = real_zero + scm_input%input_tisfc = real_zero + scm_input%input_tprcp = real_zero + scm_input%input_srflag = real_zero + scm_input%input_snwdph = real_zero + scm_input%input_shdmin = real_zero + scm_input%input_shdmax = real_zero + scm_input%input_slopetype = int_zero + scm_input%input_snoalb = real_zero + scm_input%input_sncovr = real_zero + scm_input%input_snodl = real_zero + scm_input%input_weasdl = real_zero + scm_input%input_tsfc = real_zero + scm_input%input_tsfcl = real_zero + scm_input%input_zorlw = real_zero + scm_input%input_zorll = real_zero + scm_input%input_zorli = real_zero + scm_input%input_albdirvis_lnd = real_zero + scm_input%input_albdirnir_lnd = real_zero + scm_input%input_albdifvis_lnd = real_zero + scm_input%input_albdifnir_lnd = real_zero + scm_input%input_emis_lnd = real_zero + scm_input%input_albdirvis_ice = real_zero + scm_input%input_albdifvis_ice = real_zero + scm_input%input_albdirnir_ice = real_zero + scm_input%input_albdifnir_ice = real_zero + scm_input%input_zorlwav = real_zero allocate(scm_input%input_stc(nsoil), scm_input%input_smc(nsoil), scm_input%input_slc(nsoil)) scm_input%input_stc = real_zero scm_input%input_smc = real_zero scm_input%input_slc = real_zero - allocate(scm_input%input_snicexy(nsnow),scm_input%input_snliqxy(nsnow), scm_input%input_tsnoxy(nsnow), & - scm_input%input_smoiseq(nsoil), scm_input%input_zsnsoxy(nsnow + nsoil)) - scm_input%input_snicexy = real_zero - scm_input%input_snliqxy = real_zero - scm_input%input_tsnoxy = real_zero - scm_input%input_smoiseq = real_zero - scm_input%input_zsnsoxy = real_zero - allocate(scm_input%input_tiice(nice)) scm_input%input_tiice = real_zero - allocate(scm_input%input_tslb(nsoil), scm_input%input_smois(nsoil), scm_input%input_sh2o(nsoil), & - scm_input%input_smfr(nsoil), scm_input%input_flfr(nsoil)) - scm_input%input_tslb = real_zero - scm_input%input_smois = real_zero - scm_input%input_sh2o = real_zero - scm_input%input_smfr = real_zero - scm_input%input_flfr = real_zero - - scm_input%input_tsfco = real_zero - scm_input%input_vegsrc = int_zero - scm_input%input_vegtyp = int_zero - scm_input%input_soiltyp = int_zero - scm_input%input_slopetype = int_zero - scm_input%input_vegfrac = real_zero - scm_input%input_shdmin = real_zero - scm_input%input_shdmax = real_zero - scm_input%input_zorlw = real_zero - scm_input%input_slmsk = real_zero - scm_input%input_canopy = real_zero - scm_input%input_hice = real_zero - scm_input%input_fice = real_zero - scm_input%input_tisfc = real_zero - scm_input%input_snwdph = real_zero - scm_input%input_snoalb = real_zero - scm_input%input_sncovr = real_zero - scm_input%input_area = real_zero - scm_input%input_tg3 = real_zero - scm_input%input_uustar = real_zero - scm_input%input_alvsf = real_zero - scm_input%input_alnsf = real_zero - scm_input%input_alvwf = real_zero - scm_input%input_alnwf = real_zero - scm_input%input_facsf = real_zero - scm_input%input_facwf = real_zero - scm_input%input_weasd = real_zero - scm_input%input_f10m = real_zero - scm_input%input_t2m = real_zero - scm_input%input_q2m = real_zero - scm_input%input_ffmm = real_zero - scm_input%input_ffhh = real_zero - scm_input%input_tprcp = real_zero - scm_input%input_srflag = real_zero - scm_input%input_tsfcl = real_zero - scm_input%input_zorll = real_zero - scm_input%input_zorli = real_zero - scm_input%input_zorlwav = real_zero - scm_input%input_stddev = real_zero scm_input%input_convexity = real_zero - scm_input%input_oa1 = real_zero - scm_input%input_oa2 = real_zero - scm_input%input_oa3 = real_zero - scm_input%input_oa4 = real_zero scm_input%input_ol1 = real_zero scm_input%input_ol2 = real_zero scm_input%input_ol3 = real_zero scm_input%input_ol4 = real_zero + scm_input%input_oa1 = real_zero + scm_input%input_oa2 = real_zero + scm_input%input_oa3 = real_zero + scm_input%input_oa4 = real_zero + scm_input%input_sigma = real_zero scm_input%input_theta = real_zero scm_input%input_gamma = real_zero - scm_input%input_sigma = real_zero scm_input%input_elvmax = real_zero scm_input%input_oro = real_zero scm_input%input_oro_uf = real_zero @@ -769,41 +744,44 @@ subroutine scm_input_create(scm_input, ntimes, nlev, nsoil, nsnow, nice) scm_input%input_lakefrac = real_zero scm_input%input_lakedepth = real_zero - scm_input%input_tvxy = real_zero - scm_input%input_tgxy = real_zero - scm_input%input_tahxy = real_zero - scm_input%input_canicexy = real_zero - scm_input%input_canliqxy = real_zero - scm_input%input_eahxy = real_zero - scm_input%input_cmxy = real_zero - scm_input%input_chxy = real_zero - scm_input%input_fwetxy = real_zero - scm_input%input_sneqvoxy = real_zero - scm_input%input_alboldxy = real_zero - scm_input%input_qsnowxy = real_zero - scm_input%input_wslakexy = real_zero - scm_input%input_taussxy = real_zero - scm_input%input_waxy = real_zero - scm_input%input_wtxy = real_zero - scm_input%input_zwtxy = real_zero - scm_input%input_xlaixy = real_zero - scm_input%input_xsaixy = real_zero - scm_input%input_lfmassxy = real_zero - scm_input%input_stmassxy = real_zero - scm_input%input_rtmassxy = real_zero - scm_input%input_woodxy = real_zero - scm_input%input_stblcpxy = real_zero - scm_input%input_fastcpxy = real_zero - scm_input%input_smcwtdxy = real_zero + scm_input%input_tvxy = real_zero + scm_input%input_tgxy = real_zero + scm_input%input_tahxy = real_zero + scm_input%input_canicexy = real_zero + scm_input%input_canliqxy = real_zero + scm_input%input_eahxy = real_zero + scm_input%input_cmxy = real_zero + scm_input%input_chxy = real_zero + scm_input%input_fwetxy = real_zero + scm_input%input_sneqvoxy = real_zero + scm_input%input_alboldxy = real_zero + scm_input%input_qsnowxy = real_zero + scm_input%input_wslakexy = real_zero + scm_input%input_taussxy = real_zero + scm_input%input_waxy = real_zero + scm_input%input_wtxy = real_zero + scm_input%input_zwtxy = real_zero + scm_input%input_xlaixy = real_zero + scm_input%input_xsaixy = real_zero + scm_input%input_lfmassxy = real_zero + scm_input%input_stmassxy = real_zero + scm_input%input_rtmassxy = real_zero + scm_input%input_woodxy = real_zero + scm_input%input_stblcpxy = real_zero + scm_input%input_fastcpxy = real_zero + scm_input%input_smcwtdxy = real_zero scm_input%input_deeprechxy = real_zero - scm_input%input_rechxy = real_zero - scm_input%input_snowxy = real_zero - scm_input%input_albdvis = real_zero - scm_input%input_albdnir = real_zero - scm_input%input_albivis = real_zero - scm_input%input_albinir = real_zero - scm_input%input_emiss = real_zero + scm_input%input_rechxy = real_zero + scm_input%input_snowxy = real_zero + allocate(scm_input%input_snicexy(nsnow),scm_input%input_snliqxy(nsnow), scm_input%input_tsnoxy(nsnow), & + scm_input%input_smoiseq(nsoil), scm_input%input_zsnsoxy(nsnow + nsoil)) + scm_input%input_snicexy = real_zero + scm_input%input_snliqxy = real_zero + scm_input%input_tsnoxy = real_zero + scm_input%input_smoiseq = real_zero + scm_input%input_zsnsoxy = real_zero + scm_input%input_tref = real_zero !< sea surface reference temperature for NSST (K) scm_input%input_z_c = real_zero !< sub-layer cooling thickness for NSST (m) scm_input%input_c_0 = real_zero !< coefficient 1 to calculate d(Tz)/d(Ts) for NSST @@ -833,39 +811,94 @@ subroutine scm_input_create(scm_input, ntimes, nlev, nsoil, nsnow, nice) scm_input%input_snowfallac_land = real_zero !< run-total snow accumulation on the ground over land for RUC LSM (kg m-2) scm_input%input_snowfallac_ice = real_zero !< run-total snow accumulation on the ground over ice for RUC LSM (kg m-2) scm_input%input_sncovr_ice = real_zero !< + scm_input%input_sfalb_lnd = real_zero + scm_input%input_sfalb_lnd_bck = real_zero + scm_input%input_emis_ice = real_zero scm_input%input_lai = real_zero !< leaf area index for RUC LSM - scm_input%input_sh_flux_sfc_kin = real_zero - scm_input%input_lh_flux_sfc_kin = real_zero - scm_input%input_sh_flux_sfc = real_zero - scm_input%input_lh_flux_sfc = real_zero - scm_input%input_w_ls = real_zero - scm_input%input_omega = real_zero - scm_input%input_u_g = real_zero - scm_input%input_v_g = real_zero - scm_input%input_dT_dt_rad = real_zero + allocate(scm_input%input_tslb(nsoil), scm_input%input_smois(nsoil), scm_input%input_sh2o(nsoil), & + scm_input%input_smfr(nsoil), scm_input%input_flfr(nsoil)) + scm_input%input_tslb = real_zero + scm_input%input_smois = real_zero + scm_input%input_sh2o = real_zero + scm_input%input_smfr = real_zero + scm_input%input_flfr = real_zero + + allocate(scm_input%input_pres_i(nlev+1),scm_input%input_pres_l(nlev)) + scm_input%input_pres_i = real_zero + scm_input%input_pres_l = real_zero + + allocate(scm_input%input_pres(nlev),scm_input%input_time(ntimes)) + scm_input%input_pres = real_zero + scm_input%input_time = real_zero + + allocate(scm_input%input_height(nlev), scm_input%input_thetail(nlev), scm_input%input_temp(nlev), & + scm_input%input_qt(nlev), scm_input%input_qv(nlev), scm_input%input_ql(nlev), & + scm_input%input_qi(nlev), scm_input%input_u(nlev), scm_input%input_v(nlev), scm_input%input_tke(nlev), & + scm_input%input_ozone(nlev)) + scm_input%input_height = real_zero + scm_input%input_thetail = real_zero + scm_input%input_temp = real_zero + scm_input%input_qt = real_zero + scm_input%input_qv = real_zero + scm_input%input_ql = real_zero + scm_input%input_qi = real_zero + scm_input%input_u = real_zero + scm_input%input_v = real_zero + scm_input%input_tke = real_zero + scm_input%input_ozone = real_zero + + allocate(scm_input%input_pres_surf(ntimes), scm_input%input_T_surf(ntimes)) + scm_input%input_pres_surf = real_zero + scm_input%input_T_surf = real_zero + + allocate(scm_input%input_w_ls(ntimes, nlev), scm_input%input_omega(ntimes, nlev), scm_input%input_u_g(ntimes, nlev), & + scm_input%input_v_g(ntimes, nlev), scm_input%input_u_nudge(ntimes, nlev), scm_input%input_v_nudge(ntimes, nlev), & + scm_input%input_T_nudge(ntimes, nlev), scm_input%input_thil_nudge(ntimes, nlev), scm_input%input_qt_nudge(ntimes, nlev), & + scm_input%input_dT_dt_rad(ntimes, nlev), scm_input%input_h_advec_thetail(ntimes, nlev), & + scm_input%input_h_advec_qt(ntimes, nlev), scm_input%input_v_advec_thetail(ntimes, nlev), & + scm_input%input_v_advec_qt(ntimes, nlev)) + scm_input%input_w_ls = real_zero + scm_input%input_omega = real_zero + scm_input%input_u_g = real_zero + scm_input%input_v_g = real_zero + scm_input%input_u_nudge = real_zero + scm_input%input_v_nudge = real_zero + scm_input%input_T_nudge = real_zero + scm_input%input_thil_nudge = real_zero + scm_input%input_qt_nudge = real_zero + scm_input%input_dT_dt_rad = real_zero scm_input%input_h_advec_thetail = real_zero - scm_input%input_h_advec_qt = real_zero + scm_input%input_h_advec_qt = real_zero scm_input%input_v_advec_thetail = real_zero - scm_input%input_v_advec_qt = real_zero - scm_input%input_u_nudge = real_zero - scm_input%input_v_nudge = real_zero - scm_input%input_T_nudge = real_zero - scm_input%input_thil_nudge = real_zero - scm_input%input_qt_nudge = real_zero - scm_input%input_tot_advec_t = real_zero - scm_input%input_tot_advec_theta = real_zero + scm_input%input_v_advec_qt = real_zero + + allocate(scm_input%input_sh_flux_sfc_kin(ntimes), scm_input%input_lh_flux_sfc_kin(ntimes), & + scm_input%input_sh_flux_sfc(ntimes), scm_input%input_lh_flux_sfc(ntimes)) + scm_input%input_sh_flux_sfc_kin = real_zero + scm_input%input_lh_flux_sfc_kin = real_zero + scm_input%input_sh_flux_sfc = real_zero + scm_input%input_lh_flux_sfc = real_zero + + allocate(scm_input%input_tot_advec_t(ntimes, nlev), scm_input%input_tot_advec_theta(ntimes, nlev), & + scm_input%input_tot_advec_thetal(ntimes, nlev), scm_input%input_tot_advec_qv(ntimes, nlev), & + scm_input%input_tot_advec_u(ntimes, nlev), scm_input%input_tot_advec_v(ntimes, nlev), & + scm_input%input_pres_forcing(ntimes, nlev)) + scm_input%input_tot_advec_t = real_zero + scm_input%input_tot_advec_theta = real_zero scm_input%input_tot_advec_thetal = real_zero - scm_input%input_tot_advec_qv = real_zero - scm_input%input_tot_advec_u = real_zero - scm_input%input_tot_advec_v = real_zero - scm_input%input_pres_forcing = real_zero + scm_input%input_tot_advec_qv = real_zero + scm_input%input_tot_advec_u = real_zero + scm_input%input_tot_advec_v = real_zero + scm_input%input_pres_forcing = real_zero - scm_input%input_k_T_nudge = int_one + allocate(scm_input%input_k_T_nudge(ntimes), scm_input%input_k_thil_nudge(ntimes), scm_input%input_k_qt_nudge(ntimes), & + scm_input%input_k_u_nudge(ntimes), scm_input%input_k_v_nudge(ntimes)) + scm_input%input_k_T_nudge = int_one scm_input%input_k_thil_nudge = int_one - scm_input%input_k_qt_nudge = int_one - scm_input%input_k_u_nudge = int_one - scm_input%input_k_v_nudge = int_one + scm_input%input_k_qt_nudge = int_one + scm_input%input_k_u_nudge = int_one + scm_input%input_k_v_nudge = int_one end subroutine scm_input_create @@ -988,7 +1021,7 @@ subroutine physics_set(physics, scm_input, scm_state) !this should be utilized for variables that cannot be modified or "forced" by the SCM; !most of this routine follows what is in FV3/io/FV3GFS_io.F90/sfc_prop_restart_read use scm_physical_constants, only: con_tice - use data_qc, only: conditionally_set_var + use data_qc, only: conditionally_set_var, check_missing use NetCDF_read, only: missing_value class(physics_type) :: physics @@ -998,10 +1031,38 @@ subroutine physics_set(physics, scm_input, scm_state) integer :: i,j,n real(kind=dp) :: tem1, tem logical :: missing_var(100) + real, parameter:: min_lake_orog = 200.0_dp + + !check whether input has NoahMP or RUC LSM input data + if (scm_state%model_ics .or. scm_state%lsm_ics) then + if (physics%Model%lsm == physics%Model%lsm_noahmp) then + !FV3GFS_io.F90 uses the presence of the snowxy variable in the ICs to indicate presence of NoahMP warm start + call check_missing(scm_input%input_snowxy, missing_var(1)) + if (missing_var(1)) then + physics%Model%lsm_cold_start = .true. + else + physics%Model%lsm_cold_start = .false. + end if + elseif (physics%Model%lsm == physics%Model%lsm_ruc) then + !RUC LSM uses the tslb variable as soil temperature; if it is missing, assume a cold start using Noah LSM ICs + call check_missing(scm_input%input_tslb(:), missing_var(1)) + if (missing_var(1)) then + physics%Model%lsm_cold_start = .true. + else + physics%Model%lsm_cold_start = .false. + end if + end if + end if !double check under what circumstances these should actually be set from input!!! (these overwrite the initialzation in GFS_typedefs) missing_var = .false. do i = 1, physics%Model%ncols + !since landfrac and lakefrac are read in with the orographic data (here and in FV3GFS_io), we need to set their values + !to missing when only LSM ICs are available in order for some of the logic below to work as intended. + if (scm_state%lsm_ics) then + physics%Sfcprop%landfrac(i) = missing_value + physics%Sfcprop%lakefrac(i) = missing_value + end if if (scm_state%model_ics) then write(0,'(a)') "Setting internal physics variables from the orographic section of the case input file (scalars)..." call conditionally_set_var(scm_input%input_stddev, physics%Sfcprop%hprime(i,1), "stddev", .true., missing_var(1)) @@ -1023,7 +1084,7 @@ subroutine physics_set(physics, scm_input, scm_state) call conditionally_set_var(scm_input%input_landfrac, physics%Sfcprop%landfrac(i), "landfrac", physics%Model%frac_grid, missing_var(17)) call conditionally_set_var(scm_input%input_lakefrac, physics%Sfcprop%lakefrac(i), "lakefrac", (physics%Model%lkm == 1), missing_var(18)) call conditionally_set_var(scm_input%input_lakedepth, physics%Sfcprop%lakedepth(i), "lakedepth", (physics%Model%lkm == 1), missing_var(19)) - + !write out warning if missing data for non-required variables n = 19 if ( i==1 .and. ANY( missing_var(1:n) ) ) then @@ -1040,9 +1101,9 @@ subroutine physics_set(physics, scm_input, scm_state) write(0,'(a)') "Setting internal physics variables from the surface section of the case input file (scalars)..." call conditionally_set_var(scm_input%input_slmsk, physics%Sfcprop%slmsk(i), "slmsk", (.not. physics%Model%frac_grid), missing_var(1)) call conditionally_set_var(scm_input%input_tsfco, physics%Sfcprop%tsfco(i), "tsfco", .true., missing_var(2)) - call conditionally_set_var(scm_input%input_weasd, physics%Sfcprop%weasd(i), "weasd", .false., missing_var(3)) + call conditionally_set_var(scm_input%input_weasd, physics%Sfcprop%weasd(i), "weasd", .true., missing_var(3)) call conditionally_set_var(scm_input%input_tg3, physics%Sfcprop%tg3(i), "tg3", .true., missing_var(4)) - call conditionally_set_var(scm_input%input_zorlw, physics%Sfcprop%zorlw(i), "zorlw", .true., missing_var(5)) + call conditionally_set_var(scm_input%input_zorl, physics%Sfcprop%zorl(i), "zorl", .true., missing_var(5)) call conditionally_set_var(scm_input%input_alvsf, physics%Sfcprop%alvsf(i), "alvsf", .true., missing_var(6)) call conditionally_set_var(scm_input%input_alnsf, physics%Sfcprop%alnsf(i), "alnsf", .true., missing_var(7)) call conditionally_set_var(scm_input%input_alvwf, physics%Sfcprop%alvwf(i), "alvwf", .true., missing_var(8)) @@ -1066,25 +1127,41 @@ subroutine physics_set(physics, scm_input, scm_state) call conditionally_set_var(scm_input%input_tisfc, physics%Sfcprop%tisfc(i), "tisfc", .true., missing_var(24)) call conditionally_set_var(scm_input%input_tprcp, physics%Sfcprop%tprcp(i), "tprcp", .false., missing_var(25)) call conditionally_set_var(scm_input%input_srflag, physics%Sfcprop%srflag(i), "srflag", .false., missing_var(26)) - call conditionally_set_var(scm_input%input_snwdph, physics%Sfcprop%snowd(i), "snwdph", .false., missing_var(27)) + call conditionally_set_var(scm_input%input_snwdph, physics%Sfcprop%snowd(i), "snwdph", .true., missing_var(27)) call conditionally_set_var(scm_input%input_shdmin, physics%Sfcprop%shdmin(i), "shdmin", .true., missing_var(28)) call conditionally_set_var(scm_input%input_shdmax, physics%Sfcprop%shdmax(i), "shdmax", .true., missing_var(29)) call conditionally_set_var(scm_input%input_slopetype, physics%Sfcprop%slope(i), "slopetyp", .true., missing_var(30)) call conditionally_set_var(scm_input%input_snoalb, physics%Sfcprop%snoalb(i), "snoalb", .true., missing_var(31)) call conditionally_set_var(scm_input%input_sncovr, physics%Sfcprop%sncovr(i), "sncovr", .false., missing_var(32)) - call conditionally_set_var(scm_input%input_tsfcl, physics%Sfcprop%tsfcl(i), "tsfcl", .false., missing_var(33)) - call conditionally_set_var(scm_input%input_zorll, physics%Sfcprop%zorll(i), "zorll", .false., missing_var(34)) - call conditionally_set_var(scm_input%input_zorli, physics%Sfcprop%zorli(i), "zorli", .false., missing_var(35)) + call conditionally_set_var(scm_input%input_snodl, physics%Sfcprop%snodl(i), "snodl", .false., missing_var(33)) + call conditionally_set_var(scm_input%input_weasdl, physics%Sfcprop%weasdl(i), "weasdl", .false., missing_var(34)) + call conditionally_set_var(scm_input%input_tsfc, physics%Sfcprop%tsfc(i), "tsfc", .false., missing_var(35)) + call conditionally_set_var(scm_input%input_tsfcl, physics%Sfcprop%tsfcl(i), "tsfcl", .false., missing_var(36)) + call conditionally_set_var(scm_input%input_zorlw, physics%Sfcprop%zorlw(i), "zorlw", .false., missing_var(37)) + call conditionally_set_var(scm_input%input_zorll, physics%Sfcprop%zorll(i), "zorll", .false., missing_var(38)) + call conditionally_set_var(scm_input%input_zorli, physics%Sfcprop%zorli(i), "zorli", .false., missing_var(39)) + call conditionally_set_var(scm_input%input_albdirvis_lnd, physics%Sfcprop%albdirvis_lnd(i), "albdirvis_lnd", .false., missing_var(40)) + call conditionally_set_var(scm_input%input_albdirnir_lnd, physics%Sfcprop%albdirnir_lnd(i), "albdirnir_lnd", .false., missing_var(41)) + call conditionally_set_var(scm_input%input_albdifvis_lnd, physics%Sfcprop%albdifvis_lnd(i), "albdifvis_lnd", .false., missing_var(42)) + call conditionally_set_var(scm_input%input_albdifnir_lnd, physics%Sfcprop%albdifnir_lnd(i), "albdifnir_lnd", .false., missing_var(43)) + call conditionally_set_var(scm_input%input_emis_lnd, physics%Sfcprop%emis_lnd(i), "emis_lnd", .false., missing_var(44)) + if (physics%Model%use_cice_alb .or. physics%Model%lsm == physics%Model%lsm_ruc) then + call conditionally_set_var(scm_input%input_albdirvis_ice, physics%Sfcprop%albdirvis_ice(i), "albdirvis_ice", .false., missing_var(45)) + call conditionally_set_var(scm_input%input_albdirnir_ice, physics%Sfcprop%albdirnir_ice(i), "albdirnir_ice", .false., missing_var(46)) + call conditionally_set_var(scm_input%input_albdifvis_ice, physics%Sfcprop%albdifvis_ice(i), "albdifvis_ice", .false., missing_var(47)) + call conditionally_set_var(scm_input%input_albdifnir_ice, physics%Sfcprop%albdifnir_ice(i), "albdifnir_ice", .false., missing_var(48)) + end if if (physics%Model%cplwav) then - call conditionally_set_var(scm_input%input_zorlwav, physics%Sfcprop%zorlwav(i), "zorlwav", .true., missing_var(36)) + call conditionally_set_var(scm_input%input_zorlwav, physics%Sfcprop%zorlwav(i), "zorlwav", .true., missing_var(49)) else physics%Sfcprop%zorlwav(i) = physics%Sfcprop%zorlw(i) end if + !GJF: Is this still needed? if (missing_var(32)) physics%Sfcprop%sncovr(i) = real_zero !write out warning if missing data for non-required variables - n = 36 + n = 49 if ( i==1 .and. ANY( missing_var(1:n) ) ) then write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) surface variables. This may lead to crashes or other strange behavior." write(0,'(a)') "Check scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" @@ -1097,7 +1174,7 @@ subroutine physics_set(physics, scm_input, scm_state) physics%Sfcprop%slmsk(i) = scm_state%sfc_type(i) ! tsfco is already pointing to T_surf forcing in physics_associate ! physics%Sfcprop%tsfco(i) => scm_state%T_surf - physics%Sfcprop%zorlw(i) = scm_state%sfc_roughness_length_cm(i) + physics%Sfcprop%zorl(i) = scm_state%sfc_roughness_length_cm(i) ! tisfc is already pointing to T_surf forcing in physics_associate ! physics%Sfcprop%tisfc(i) => scm_state%T_surf ! tsfcl is already pointing to T_surf forcing in physics_associate @@ -1110,62 +1187,104 @@ subroutine physics_set(physics, scm_input, scm_state) ! physics%Model%ivegsrc = scm_input%input_vegsrc !end if - if(scm_state%model_ics .and. physics%Model%frac_grid) then ! obtain slmsk from landfrac - !! next 5 lines are temporary till lake model is available - !if (physics%Sfcprop%lakefrac(i) > real_zero) then - ! !physics%Sfcprop%lakefrac(i) = nint(physics%Sfcprop%lakefrac(i)) - ! physics%Sfcprop%landfrac(i) = real_one - physics%Sfcprop%lakefrac(i) - ! if (physics%Sfcprop%lakefrac(i) == real_zero) physics%Sfcprop%fice(i) = real_zero - !endif - !physics%Sfcprop%slmsk(i) = ceiling(physics%Sfcprop%landfrac(i)) - !if (physics%Sfcprop%fice(i) > physics%Model%min_lakeice .and. physics%Sfcprop%landfrac(i) == real_zero) physics%Sfcprop%slmsk(i) = 2 ! land dominates ice if co-exist - physics%Sfcprop%slmsk(i) = ceiling(physics%Sfcprop%landfrac(i)) !nint/floor are options - else ! obtain landfrac from slmsk - if (physics%Sfcprop%slmsk(i) > 1.9_dp) then + if(scm_state%model_ics .or. scm_state%lsm_ics) then + if (physics%Sfcprop%stype(i) == 14 .or. physics%Sfcprop%stype(i)+0.5 <= 0) then physics%Sfcprop%landfrac(i) = real_zero - else - physics%Sfcprop%landfrac(i) = physics%Sfcprop%slmsk(i) + physics%Sfcprop%stype(i) = 0 + if (physics%Sfcprop%lakefrac(i) > real_zero) then + physics%Sfcprop%lakefrac(i) = real_one + endif endif - endif - - if (physics%Sfcprop%lakefrac(i) > real_zero) then - physics%Sfcprop%oceanfrac(i) = real_zero ! lake & ocean don't coexist in a cell - !GJF: the following code was taken from fv3atm/develop branch commit d10e450 from 2020/12/2, but the if statements seem like they are in error (reproduced here) - if (physics%Sfcprop%slmsk(i) /= real_one) then - if (physics%Sfcprop%fice(i) >= physics%Model%min_lakeice) then - if (physics%Sfcprop%slmsk(i) < 1.9_dp) & - write(*,'(a,2i3,3f6.2)') 'reset lake slmsk=2 at i=' & - ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%lakefrac(i) - physics%Sfcprop%slmsk(i) = 2. - else if (physics%Sfcprop%slmsk(i) > 1.e-7) then - write(*,'(a,2i3,3f6.2)') 'reset lake slmsk=0 at i=' & - ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%lakefrac(i) - physics%Sfcprop%slmsk(i) = real_zero - end if - end if - !if (physics%Sfcprop%fice(i) < physics%Model%min_lakeice) then - ! physics%Sfcprop%fice(i) = real_zero - ! if (physics%Sfcprop%slmsk(i) == 2) physics%Sfcprop%slmsk(i) = 0 - !endif - else - physics%Sfcprop%oceanfrac(i) = real_one - physics%Sfcprop%landfrac(i) - if (physics%Sfcprop%slmsk(i) /= real_one) then - if (physics%Sfcprop%fice(i) >= physics%Model%min_seaice) then - if (physics%Sfcprop%slmsk(i) < 1.9_dp) & - write(*,'(a,2i3,3f6.2)') 'reset sea slmsk=2 at i,=' & - ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%landfrac(i) - physics%Sfcprop%slmsk(i) = 2. - else if (physics%Sfcprop%slmsk(i) > 1.e-7) then - write(*,'(a,2i3,4f6.2)') 'reset sea slmsk=0 at i=' & - ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%landfrac(i) - physics%Sfcprop%slmsk(i) = real_zero - end if - end if - ! if (physics%Sfcprop%fice(i) < physics%Model%min_seaice) then - ! physics%Sfcprop%fice(i) = real_zero - ! if (physics%Sfcprop%slmsk(i) == 2) physics%Sfcprop%slmsk(i) = 0 - ! endif - endif + + if (physics%Model%frac_grid) then + if (physics%Sfcprop%landfrac(i) > -999.0_dp) then + physics%Sfcprop%slmsk(i) = ceiling(physics%Sfcprop%landfrac(i)-1.0e-6) + if (physics%Sfcprop%slmsk(i) == 1 .and. physics%Sfcprop%stype(i) == 14) & + physics%Sfcprop%slmsk(i) = 0 + if (physics%Sfcprop%lakefrac(i) > real_zero) then + physics%Sfcprop%oceanfrac(i) = real_zero ! lake & ocean don't coexist in a cell + if (nint(physics%Sfcprop%slmsk(i)) /= 1) then + if(physics%Sfcprop%fice(i) >= physics%Model%min_lakeice) then + physics%Sfcprop%slmsk(i) = 2 + else + physics%Sfcprop%slmsk(i) = 0 + endif + endif + else + physics%Sfcprop%lakefrac(i) = real_zero + physics%Sfcprop%oceanfrac(i) = real_one - physics%Sfcprop%landfrac(i) + if (nint(physics%Sfcprop%slmsk(i)) /= 1) then + if (physics%Sfcprop%fice(i) >= physics%Model%min_seaice) then + physics%Sfcprop%slmsk(i) = 2 + else + physics%Sfcprop%slmsk(i) = 0 + endif + endif + endif + else + physics%Model%frac_grid = .false. + if (nint(physics%Sfcprop%slmsk(i)) == 1) then + physics%Sfcprop%landfrac(i) = real_one + physics%Sfcprop%lakefrac(i) = real_zero + physics%Sfcprop%oceanfrac(i) = real_zero + else + if (physics%Sfcprop%slmsk(i) < 0.1_dp .or. physics%Sfcprop%slmsk(i) > 1.9_dp) then + physics%Sfcprop%landfrac(i) = real_zero + if (physics%Sfcprop%oro_uf(i) > min_lake_orog) then ! lakes + physics%Sfcprop%lakefrac(i) = real_one + physics%Sfcprop%oceanfrac(i) = real_zero + else ! ocean + physics%Sfcprop%lakefrac(i) = real_zero + physics%Sfcprop%oceanfrac(i) = real_one + endif + endif + endif + endif + else ! not a fractional grid + if (physics%Sfcprop%landfrac(i) > -999.0_dp) then + if (physics%Sfcprop%lakefrac(i) > real_zero) then + physics%Sfcprop%oceanfrac(i) = real_zero + physics%Sfcprop%landfrac(i) = real_zero + physics%Sfcprop%lakefrac(i) = real_one + physics%Sfcprop%slmsk(i) = real_zero + if (physics%Sfcprop%fice(i) >= physics%Model%min_lakeice) physics%Sfcprop%slmsk(i) = 2.0 + else + physics%Sfcprop%slmsk(i) = nint(physics%Sfcprop%landfrac(i)) + if (physics%Sfcprop%stype(i) <= 0 .or. physics%Sfcprop%stype(i) == 14) & + physics%Sfcprop%slmsk(i) = real_zero + if (nint(physics%Sfcprop%slmsk(i)) == 0) then + physics%Sfcprop%oceanfrac(i) = real_one + physics%Sfcprop%landfrac(i) = real_zero + physics%Sfcprop%lakefrac(i) = real_zero + if (physics%Sfcprop%fice(i) >= physics%Model%min_seaice) physics%Sfcprop%slmsk(i) = 2.0 + else + physics%Sfcprop%landfrac(i) = real_one + physics%Sfcprop%lakefrac(i) = real_zero + physics%Sfcprop%oceanfrac(i) = real_zero + endif + endif + else + if (nint(physics%Sfcprop%slmsk(i)) == 1 .and. physics%Sfcprop%stype(i) > 0 & + .and. physics%Sfcprop%stype(i) /= 14) then + physics%Sfcprop%landfrac(i) = real_one + physics%Sfcprop%lakefrac(i) = real_zero + physics%Sfcprop%oceanfrac(i) = real_zero + else + physics%Sfcprop%slmsk(i) = real_zero + physics%Sfcprop%landfrac(i) = real_zero + if (physics%Sfcprop%oro_uf(i) > min_lake_orog) then ! lakes + physics%Sfcprop%lakefrac(i) = real_one + physics%Sfcprop%oceanfrac(i) = real_zero + if (physics%Sfcprop%fice(i) > physics%Model%min_lakeice) physics%Sfcprop%slmsk(i) = 2.0 + else ! ocean + physics%Sfcprop%lakefrac(i) = real_zero + physics%Sfcprop%oceanfrac(i) = real_one + if (physics%Sfcprop%fice(i) > physics%Model%min_seaice) physics%Sfcprop%slmsk(i) = 2.0 + endif + endif + endif + endif + end if !--- NSSTM variables if (physics%Model%nstf_name(1) > 0) then @@ -1214,31 +1333,35 @@ subroutine physics_set(physics, scm_input, scm_state) endif endif - if ((scm_state%model_ics .or. scm_state%lsm_ics) .and. physics%Model%lsm == physics%Model%lsm_ruc) then !.and. warm_start (not implemented here -- assuming that RUC LSM has warm start data from file) + if ((scm_state%model_ics .or. scm_state%lsm_ics) .and. physics%Model%lsm == physics%Model%lsm_ruc .and. .not. physics%Model%lsm_cold_start) then !--- Extra RUC LSM variables write(0,'(a)') "Setting internal physics variables from the RUC LSM section of the case input file (scalars)..." - call conditionally_set_var(scm_input%input_wetness, physics%Sfcprop%wetness(i), "wetness", .false., missing_var(1)) - call conditionally_set_var(scm_input%input_clw_surf_land, physics%Sfcprop%clw_surf_land(i), "clw_surf_land", .false., missing_var(2)) - call conditionally_set_var(scm_input%input_clw_surf_ice, physics%Sfcprop%clw_surf_ice(i), "clw_surf_ice", .false., missing_var(3)) - call conditionally_set_var(scm_input%input_qwv_surf_land, physics%Sfcprop%qwv_surf_land(i), "qwv_surf_land", .false., missing_var(4)) - call conditionally_set_var(scm_input%input_qwv_surf_ice, physics%Sfcprop%qwv_surf_ice(i), "qwv_surf_ice", .false., missing_var(5)) - call conditionally_set_var(scm_input%input_tsnow_land, physics%Sfcprop%tsnow_land(i), "tsnow_land", .false., missing_var(6)) - call conditionally_set_var(scm_input%input_tsnow_ice, physics%Sfcprop%tsnow_ice(i), "tsnow_ice", .false., missing_var(7)) - call conditionally_set_var(scm_input%input_snowfallac_land, physics%Sfcprop%snowfallac_land(i), "snowfallac_land", .false., missing_var(8)) - call conditionally_set_var(scm_input%input_snowfallac_ice, physics%Sfcprop%snowfallac_ice(i), "snowfallac_ice", .false., missing_var(9)) + call conditionally_set_var(scm_input%input_wetness, physics%Sfcprop%wetness(i), "wetness", .true., missing_var(1)) + call conditionally_set_var(scm_input%input_clw_surf_land, physics%Sfcprop%clw_surf_land(i), "clw_surf_land", .true., missing_var(2)) + call conditionally_set_var(scm_input%input_clw_surf_ice, physics%Sfcprop%clw_surf_ice(i), "clw_surf_ice", .true., missing_var(3)) + call conditionally_set_var(scm_input%input_qwv_surf_land, physics%Sfcprop%qwv_surf_land(i), "qwv_surf_land", .true., missing_var(4)) + call conditionally_set_var(scm_input%input_qwv_surf_ice, physics%Sfcprop%qwv_surf_ice(i), "qwv_surf_ice", .true., missing_var(5)) + call conditionally_set_var(scm_input%input_tsnow_land, physics%Sfcprop%tsnow_land(i), "tsnow_land", .true., missing_var(6)) + call conditionally_set_var(scm_input%input_tsnow_ice, physics%Sfcprop%tsnow_ice(i), "tsnow_ice", .true., missing_var(7)) + call conditionally_set_var(scm_input%input_snowfallac_land, physics%Sfcprop%snowfallac_land(i), "snowfallac_land", .true., missing_var(8)) + call conditionally_set_var(scm_input%input_snowfallac_ice, physics%Sfcprop%snowfallac_ice(i), "snowfallac_ice", .true., missing_var(9)) call conditionally_set_var(scm_input%input_sncovr_ice, physics%Sfcprop%sncovr_ice(i), "sncovr_ice", .false., missing_var(10)) + call conditionally_set_var(scm_input%input_sfalb_lnd, physics%Sfcprop%sfalb_lnd(i), "sfalb_lnd", .true., missing_var(11)) + call conditionally_set_var(scm_input%input_sfalb_lnd_bck, physics%Sfcprop%sfalb_lnd_bck(i), "sfalb_lnd_bck", .true., missing_var(12)) + call conditionally_set_var(scm_input%input_sfalb_ice, physics%Sfcprop%sfalb_ice(i), "sfalb_ice", .true., missing_var(13)) + call conditionally_set_var(scm_input%input_emis_ice, physics%Sfcprop%emis_ice(i), "emis_ice", .true., missing_var(14)) if (physics%Model%lsm == physics%Model%lsm_ruc .and. physics%Model%rdlai) then !when rdlai = T, RUC LSM expects the LAI to be read in, hence the required variable attribute below - call conditionally_set_var(scm_input%input_lai, physics%Sfcprop%xlaixy(i), "lai", .true., missing_var(11)) + call conditionally_set_var(scm_input%input_lai, physics%Sfcprop%xlaixy(i), "lai", .true., missing_var(15)) end if !if sncovr_ice is missing, set to the land value if(missing_var(10)) then - call conditionally_set_var(scm_input%input_sncovr, physics%Sfcprop%sncovr_ice(i), "sncovr_ice", .false., missing_var(10)) + call conditionally_set_var(scm_input%input_sncovr, physics%Sfcprop%sncovr_ice(i), "sncovr_ice", .true., missing_var(10)) end if !write out warning if missing data for non-required variables - n = 11 + n = 15 if ( i==1 .and. ANY( missing_var(1:n) ) ) then write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) surface variables for RUC LSM. This may lead to crashes or other strange behavior." write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" @@ -1247,6 +1370,23 @@ subroutine physics_set(physics, scm_input, scm_state) end do end if missing_var = .false. + elseif ((scm_state%model_ics .or. scm_state%lsm_ics) .and. physics%Model%lsm == physics%Model%lsm_ruc .and. physics%Model%lsm_cold_start) then + call conditionally_set_var(scm_input%input_sncovr, physics%Sfcprop%sncovr_ice(i), "sncovr_ice", .true., missing_var(1)) + if (physics%Model%rdlai) then + !when rdlai = T, RUC LSM expects the LAI to be read in, hence the required variable attribute below + call conditionally_set_var(scm_input%input_lai, physics%Sfcprop%xlaixy(i), "lai", .true., missing_var(2)) + end if + + !write out warning if missing data for non-required variables + n = 2 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) surface variables for RUC LSM. This may lead to crashes or other strange behavior." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. elseif ((scm_state%model_ics .or. scm_state%lsm_ics) .and. physics%Model%lsm == physics%Model%lsm_noahmp) then write(0,'(a)') "Setting internal physics variables from the NoahMP section of the case input file (scalars)..." !all of these can be missing, since a method exists to "cold start" these variables @@ -1279,14 +1419,9 @@ subroutine physics_set(physics, scm_input, scm_state) call conditionally_set_var(scm_input%input_smcwtdxy, physics%Sfcprop%smcwtdxy(i), "smcwtdxy", .false., missing_var(27)) call conditionally_set_var(scm_input%input_deeprechxy, physics%Sfcprop%deeprechxy(i), "deeprechxy", .false., missing_var(28)) call conditionally_set_var(scm_input%input_rechxy, physics%Sfcprop%rechxy(i), "rechxy", .false., missing_var(29)) - call conditionally_set_var(scm_input%input_albdvis, physics%Sfcprop%albdirvis_lnd(i), "albdirvis_lnd", .false., missing_var(30)) - call conditionally_set_var(scm_input%input_albdnir, physics%Sfcprop%albdirnir_lnd(i), "albdirnir_lnd", .false., missing_var(31)) - call conditionally_set_var(scm_input%input_albivis, physics%Sfcprop%albdifvis_lnd(i), "albdifvis_lnd", .false., missing_var(32)) - call conditionally_set_var(scm_input%input_albinir, physics%Sfcprop%albdirnir_lnd(i), "albdifnir_lnd", .false., missing_var(33)) - call conditionally_set_var(scm_input%input_emiss, physics%Sfcprop%emis_lnd(i), "emis_lnd", .false., missing_var(34)) !write out warning if missing data for non-required variables - n = 34 + n = 29 if ( i==1 .and. ANY( missing_var(1:n) ) ) then write(0,'(a)') "INPUT CHECK: Some missing input data was found related to surface variables for NoahMP LSM. Due to this, a cold-start algorithm to initialize variables will be used." write(0,'(a)') "Check scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" @@ -1298,36 +1433,19 @@ subroutine physics_set(physics, scm_input, scm_state) end if if ((scm_state%model_ics .or. scm_state%lsm_ics) .and. (physics%Model%lsm == physics%Model%lsm_noah .or. & - physics%Model%lsm == physics%Model%lsm_noahmp)) then !.or. (.not. warm_start) from FV3GFS_io is not implemented - !check for nonmissing values - !--- 3D variables - ! do lsoil = 1,physics%Model%lsoil - ! physics%Sfcprop%stc(i,lsoil) = scm_input%input_stc(lsoil) !--- stc - ! physics%Sfcprop%smc(i,lsoil) = scm_input%input_smc(lsoil) !--- smc - ! physics%Sfcprop%slc(i,lsoil) = scm_input%input_slc(lsoil) !--- slc - ! enddo + physics%Model%lsm == physics%Model%lsm_noahmp .or. physics%Model%lsm_cold_start)) then + call conditionally_set_var(scm_input%input_stc(:), physics%Sfcprop%stc(i,:), "stc", .true., missing_var(1)) call conditionally_set_var(scm_input%input_smc(:), physics%Sfcprop%smc(i,:), "smc", .true., missing_var(2)) call conditionally_set_var(scm_input%input_slc(:), physics%Sfcprop%slc(i,:), "slc", .true., missing_var(3)) if (physics%Model%lsm == physics%Model%lsm_noahmp) then - ! do lsoil = -2, 0 - ! physics%Sfcprop%snicexy(i,lsoil) = scm_input%input_snicexy(lsoil) - ! physics%Sfcprop%snliqxy(i,lsoil) = scm_input%input_snliqxy(lsoil) - ! physics%Sfcprop%tsnoxy(i,lsoil) = scm_input%input_tsnoxy(lsoil) - ! enddo call conditionally_set_var(scm_input%input_snicexy(:), physics%Sfcprop%snicexy(i,:), "snicexy", .false., missing_var(4)) call conditionally_set_var(scm_input%input_snliqxy(:), physics%Sfcprop%snliqxy(i,:), "sliqexy", .false., missing_var(5)) call conditionally_set_var(scm_input%input_tsnoxy(:), physics%Sfcprop%tsnoxy(i,:), "tsnoxy", .false., missing_var(6)) - ! do lsoil = 1, 4 - ! physics%Sfcprop%smoiseq(i,lsoil) = scm_input%input_smoiseq(lsoil) - ! enddo call conditionally_set_var(scm_input%input_smoiseq(:), physics%Sfcprop%smoiseq(i,:), "smoiseq", .false., missing_var(7)) - ! do lsoil = -2, 4 - ! physics%Sfcprop%zsnsoxy(i,lsoil) = scm_input%input_zsnsoxy(lsoil) - ! enddo call conditionally_set_var(scm_input%input_zsnsoxy(:), physics%Sfcprop%zsnsoxy(i,:), "zsnzoxy", .false., missing_var(8)) !write out warning if missing data for non-required variables @@ -1343,14 +1461,6 @@ subroutine physics_set(physics, scm_input, scm_state) endif else if ((scm_state%model_ics .or. scm_state%lsm_ics) .and. physics%Model%lsm == physics%Model%lsm_ruc) then - !--- 3D variables - ! do lsoil = 1,Model%lsoil_lsm - ! physics%Sfcprop%tslb(i,lsoil) = scm_input%input_tslb(lsoil) - ! physics%Sfcprop%smois(i,lsoil) = scm_input%input_smois(lsoil) - ! physics%Sfcprop%sh2o(i,lsoil) = scm_input%input_sh2o(lsoil) - ! physics%Sfcprop%keepsmfr(i,lsoil) = scm_input%input_smfr(lsoil) - ! physics%Sfcprop%flag_frsoil(i,lsoil) = scm_input%input_flfr(lsoil) - ! enddo call conditionally_set_var(scm_input%input_tslb(:), physics%Sfcprop%tslb(i,:), "tslb", .false., missing_var(1)) call conditionally_set_var(scm_input%input_smois(:), physics%Sfcprop%smois(i,:), "smois", .false., missing_var(2)) call conditionally_set_var(scm_input%input_sh2o(:), physics%Sfcprop%sh2o(i,:), "sh2o", .false., missing_var(3)) @@ -1371,67 +1481,102 @@ subroutine physics_set(physics, scm_input, scm_state) if (scm_state%model_ics .or. scm_state%lsm_ics) then !check for nonmissing values - ! do k = 1,Model%kice - ! physics%Sfcprop%tiice(i,k) = scm_input%input_tiice(k) !--- internal ice temp - ! enddo call conditionally_set_var(scm_input%input_tiice(:), physics%Sfcprop%tiice(i,:), "tiice", .false., missing_var(1)) if (missing_var(1)) then write(0,'(a)') "INPUT CHECK: Some missing input data was found related to the internal sea ice temperature. These will be set from the internal soil temperature variable (stc)." end if end if + if (scm_state%model_ics .or. scm_state%lsm_ics) then + if (scm_input%input_snodl <= real_zero) then + if (physics%Sfcprop%landfrac(i) > real_zero) then + tem = real_one / physics%Sfcprop%landfrac(i) + physics%Sfcprop%snodl(i) = physics%Sfcprop%snowd(i) * tem + else + physics%Sfcprop%snodl(i) = real_zero + endif + end if + + if (scm_input%input_weasdl <= real_zero) then + if (physics%Sfcprop%landfrac(i) > real_zero) then + tem = real_one / physics%Sfcprop%landfrac(i) + physics%Sfcprop%weasdl(i) = physics%Sfcprop%weasd(i) * tem + else + physics%Sfcprop%weasdl(i) = real_zero + endif + end if + end if + if (scm_input%input_tsfcl <= real_zero) then physics%Sfcprop%tsfcl(i) = physics%Sfcprop%tsfco(i) !--- compute tsfcl from existing variables end if + if (scm_input%input_zorlw <= real_zero) then + physics%Sfcprop%zorlw(i) = physics%Sfcprop%zorl(i) !--- compute zorlw from existing variables + end if + if (scm_input%input_zorll <= real_zero) then - physics%Sfcprop%zorll(i) = physics%Sfcprop%zorlw(i) !--- compute zorll from existing variables + physics%Sfcprop%zorll(i) = physics%Sfcprop%zorl(i) !--- compute zorll from existing variables end if if (scm_input%input_zorli <= real_zero) then - physics%Sfcprop%zorli(i) = physics%Sfcprop%zorlw(i) !--- compute zorli from existing variables + physics%Sfcprop%zorli(i) = physics%Sfcprop%zorl(i) !--- compute zorli from existing variables end if + if (physics%Model%use_cice_alb) then + if (scm_input%input_albdirvis_ice <= real_zero) then + if (physics%Sfcprop%oceanfrac(i) > real_zero .and. & + physics%Sfcprop%fice(i) >= physics%Model%min_seaice) then + physics%Sfcprop%albdirvis_ice(i) = 0.6_dp + physics%Sfcprop%albdifvis_ice(i) = 0.6_dp + physics%Sfcprop%albdirnir_ice(i) = 0.6_dp + physics%Sfcprop%albdifnir_ice(i) = 0.6_dp + endif + endif + endif + if (scm_input%input_zorlwav <= real_zero) then physics%Sfcprop%zorlwav(i) = physics%Sfcprop%zorlw(i) !--- compute zorlwav from existing variables end if - if(physics%Model%frac_grid .and. (scm_state%model_ics .or. scm_state%lsm_ics)) then ! 3-way composite - if( physics%Model%phour < 1.e-7) physics%Sfcprop%tsfco(i) = max(con_tice, physics%Sfcprop%tsfco(i)) - tem1 = real_one - physics%Sfcprop%landfrac(i) - tem = tem1 * physics%Sfcprop%fice(i) ! tem = ice fraction wrt whole cell - physics%Sfcprop%zorl(i) = physics%Sfcprop%zorll(i) * physics%Sfcprop%landfrac(i) & - + physics%Sfcprop%zorli(i) * tem & - + physics%Sfcprop%zorlw(i) * (tem1-tem) - - physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfcl(i) * physics%Sfcprop%landfrac(i) & - + physics%Sfcprop%tisfc(i) * tem & - + physics%Sfcprop%tsfco(i) * (tem1-tem) - else - !--- specify tsfcl/zorll/zorli from existing variable tsfco/zorlw - ! physics%Sfcprop%tsfcl(i) = physics%Sfcprop%tsfco(i) - ! physics%Sfcprop%zorll(i) = physics%Sfcprop%zorlw(i) - ! physics%Sfcprop%zorli(i) = physics%Sfcprop%zorlw(i) - ! physics%Sfcprop%zorl(i) = physics%Sfcprop%zorlw(i) - ! physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfco(i) - if (physics%Sfcprop%slmsk(i) == 1) then - physics%Sfcprop%zorl(i) = physics%Sfcprop%zorll(i) - physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfcl(i) - else - tem = real_one - physics%Sfcprop%fice(i) - physics%Sfcprop%zorl(i) = physics%Sfcprop%zorli(i) * physics%Sfcprop%fice(i) & - + physics%Sfcprop%zorlw(i) * tem - - physics%Sfcprop%tsfc(i) = physics%Sfcprop%tisfc(i) * physics%Sfcprop%fice(i) & - + physics%Sfcprop%tsfco(i) * tem - endif - endif ! if (Model%frac_grid) + if (physics%Model%lsm_cold_start) then + if(physics%Model%frac_grid .and. (scm_state%model_ics .or. scm_state%lsm_ics)) then ! 3-way composite + if( physics%Model%phour < 1.e-7) physics%Sfcprop%tsfco(i) = max(con_tice, physics%Sfcprop%tsfco(i)) + tem1 = real_one - physics%Sfcprop%landfrac(i) + tem = tem1 * physics%Sfcprop%fice(i) ! tem = ice fraction wrt whole cell + physics%Sfcprop%zorl(i) = physics%Sfcprop%zorll(i) * physics%Sfcprop%landfrac(i) & + + physics%Sfcprop%zorli(i) * tem & + + physics%Sfcprop%zorlw(i) * (tem1-tem) + + physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfcl(i) * physics%Sfcprop%landfrac(i) & + + physics%Sfcprop%tisfc(i) * tem & + + physics%Sfcprop%tsfco(i) * (tem1-tem) + else + !--- specify tsfcl/zorll/zorli from existing variable tsfco/zorlw + ! physics%Sfcprop%tsfcl(i) = physics%Sfcprop%tsfco(i) + ! physics%Sfcprop%zorll(i) = physics%Sfcprop%zorlw(i) + ! physics%Sfcprop%zorli(i) = physics%Sfcprop%zorlw(i) + ! physics%Sfcprop%zorl(i) = physics%Sfcprop%zorlw(i) + ! physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfco(i) + if (physics%Sfcprop%slmsk(i) == 1) then + physics%Sfcprop%zorl(i) = physics%Sfcprop%zorll(i) + physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfcl(i) + else + tem = real_one - physics%Sfcprop%fice(i) + physics%Sfcprop%zorl(i) = physics%Sfcprop%zorli(i) * physics%Sfcprop%fice(i) & + + physics%Sfcprop%zorlw(i) * tem + + physics%Sfcprop%tsfc(i) = physics%Sfcprop%tisfc(i) * physics%Sfcprop%fice(i) & + + physics%Sfcprop%tsfco(i) * tem + endif + endif ! if (Model%frac_grid) + endif !if (physics%Model%lsm_cold_start) if ((scm_state%model_ics .or. scm_state%lsm_ics) .and. MAXVAL(scm_input%input_tiice) < real_zero) then physics%Sfcprop%tiice(i,1) = physics%Sfcprop%stc(i,1) !--- initialize internal ice temp from soil temp at layer 1 physics%Sfcprop%tiice(i,2) = physics%Sfcprop%stc(i,2) !--- initialize internal ice temp from soil temp at layer 2 end if - + end do end subroutine physics_set