diff --git a/physics/GFS_calpreciptype.f90 b/physics/GFS_calpreciptype.f90 index 3fe499c9b..eb3951391 100644 --- a/physics/GFS_calpreciptype.f90 +++ b/physics/GFS_calpreciptype.f90 @@ -75,7 +75,9 @@ end subroutine GFS_calpreciptype_finalize !! | snow0 | lwe_thickness_of_snow_amount_per_day | snow fall over 24h period | mm | 1 | real | kind_phys | in | F | !! | graupel0 | lwe_thickness_of_graupel_amount_per_day | graupel fall over 24h period | mm | 1 | real | kind_phys | in | F | !! | cplflx | flag_for_flux_coupling | flag controlling cplflx collection (default off) | flag | 0 | logical | | in | F | +!! | cplchm | flag_for_chemistry_coupling | flag controlling cplchm collection (default off) | flag | 0 | logical | | in | F | !! | rain_cpl | lwe_thickness_of_precipitation_amount_for_coupling | total rain precipitation for model coupling | m | 1 | real | kind_phys | inout | F | +!! | rainc_cpl | lwe_thickness_of_convective_precipitation_amount_for_coupling | total convective precipitation | m | 1 | real | kind_phys | inout | F | !! | snow_cpl | lwe_thickness_of_snow_amount_for_coupling | total snow precipitation for model coupling | m | 1 | real | kind_phys | inout | F | !! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | @@ -91,7 +93,8 @@ subroutine GFS_calpreciptype_run(kdt,nrcm,im,ix,lm,lp1,randomno, & totprcpb,toticeb,totsnwb,totgrpb, & hrate,mrate,save_t,save_qv, & rain0,ice0,snow0,graupel0, & - cplflx,rain_cpl,snow_cpl,errmsg,errflg) + cplflx,cplchm,rain_cpl,rainc_cpl, & + snow_cpl,errmsg,errflg) !$$$ subprogram documentation block ! . . . @@ -134,9 +137,9 @@ subroutine GFS_calpreciptype_run(kdt,nrcm,im,ix,lm,lp1,randomno, & real(kind=kind_phys), dimension(ix,lm), intent(in) :: save_t,save_qv ! rain0, ice0, snow0, graupel0 only allocated if GFDL MP is selected (imp_physics == imp_physics_gfdl) real(kind=kind_phys), dimension(:), intent(in) :: rain0,ice0,snow0,graupel0 - logical, intent(in) :: cplflx + logical, intent(in) :: cplflx,cplchm ! rain_cpl, snow_cpl only allocated if cplflx == .true. - real(kind=kind_phys), dimension(:), intent(inout) :: rain_cpl,snow_cpl + real(kind=kind_phys), dimension(:), intent(inout) :: rain_cpl,rainc_cpl,snow_cpl character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -399,7 +402,8 @@ subroutine GFS_calpreciptype_run(kdt,nrcm,im,ix,lm,lp1,randomno, & crain = 0.0 csnow = rainc(i) endif - if ((snow0(i)+ice0(i)+graupel0(i)+csnow) > (rain0(i)+crain)) then +! if ((snow0(i)+ice0(i)+graupel0(i)+csnow) > (rain0(i)+crain)) then + if ((snow0(i)+ice0(i)+graupel0(i)+csnow) > 0.0) then srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) endif enddo @@ -425,6 +429,13 @@ subroutine GFS_calpreciptype_run(kdt,nrcm,im,ix,lm,lp1,randomno, & enddo endif + if ((cplchm).and.(.not.cplflx)) then + do i = 1, im + rain_cpl(i) = rain_cpl(i) + rain(i) + rainc_cpl(i) = rainc_cpl(i) + rainc(i) + enddo + endif + return end subroutine GFS_calpreciptype_run ! diff --git a/physics/GFS_phys_time_vary.fv3.f90 b/physics/GFS_phys_time_vary.fv3.f90 index 3355123ff..aa38ed4ef 100644 --- a/physics/GFS_phys_time_vary.fv3.f90 +++ b/physics/GFS_phys_time_vary.fv3.f90 @@ -1,20 +1,174 @@ !> \file GFS_phys_time_vary.f90 !! Contains code related to GFS physics suite setup (physics part of time_vary_step) - module GFS_phys_time_vary + module GFS_phys_time_vary + + use ozne_def, only : levozp, oz_coeff, oz_lat, oz_pres, oz_time, ozplin + use ozinterp, only : read_o3data, setindxoz, ozinterpol + + use h2o_def, only : levh2o, h2o_coeff, h2o_lat, h2o_pres, h2o_time, h2oplin + use h2ointerp, only : read_h2odata, setindxh2o, h2ointerpol + + implicit none private public GFS_phys_time_vary_init, GFS_phys_time_vary_run, GFS_phys_time_vary_finalize + logical :: is_initialized = .false. + contains - subroutine GFS_phys_time_vary_init () +!> \section arg_table_GFS_phys_time_vary_init Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | Data | FV3-GFS_Data_type_all_blocks | Fortran DDT containing FV3-GFS data | DDT | 1 | GFS_data_type | | inout | F | +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_control_type | | inout | F | +!! | Interstitial | FV3-GFS_Interstitial_type_all_threads | Fortran DDT containing FV3-GFS interstitial data | DDT | 1 | GFS_interstitial_type | | inout | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! + subroutine GFS_phys_time_vary_init (Data, Model, Interstitial, errmsg, errflg) + + use GFS_typedefs, only: GFS_control_type, GFS_data_type, GFS_interstitial_type + + implicit none + + ! Interface variables + type(GFS_data_type), intent(inout) :: Data(:) + type(GFS_control_type), intent(inout) :: Model + type(GFS_interstitial_type), intent(inout) :: Interstitial(:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: nb, nblks, nt, nthrds + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (is_initialized) return + + nblks = size(Model%blksz) + nthrds = size(Interstitial) + + call read_o3data (Model%ntoz, Model%me, Model%master) + + ! Consistency check that the hardcoded values for levozp and + ! oz_coeff in GFS_typedefs.F90 match what is set by read_o3data + ! in GFS_typedefs.F90: allocate (Tbd%ozpl (IM,levozp,oz_coeff)) + if (size(Data(1)%Tbd%ozpl, dim=2).ne.levozp) then + write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ", & + "levozp from read_o3data does not match value in GFS_typedefs.F90: ", & + levozp, " /= ", size(Data(1)%Tbd%ozpl, dim=2) + errflg = 1 + return + end if + if (size(Data(1)%Tbd%ozpl, dim=3).ne.oz_coeff) then + write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ", & + "oz_coeff from read_o3data does not match value in GFS_typedefs.F90: ", & + oz_coeff, " /= ", size(Data(1)%Tbd%ozpl, dim=3) + errflg = 1 + return + end if + + call read_h2odata (Model%h2o_phys, Model%me, Model%master) + + ! Consistency check that the hardcoded values for levh2o and + ! h2o_coeff in GFS_typedefs.F90 match what is set by read_o3data + ! in GFS_typedefs.F90: allocate (Tbd%h2opl (IM,levh2o,h2o_coeff)) + if (size(Data(1)%Tbd%h2opl, dim=2).ne.levh2o) then + write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ", & + "levh2o from read_h2odata does not match value in GFS_typedefs.F90: ", & + levh2o, " /= ", size(Data(1)%Tbd%h2opl, dim=2) + errflg = 1 + return + end if + if (size(Data(1)%Tbd%h2opl, dim=3).ne.h2o_coeff) then + write(errmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ", & + "h2o_coeff from read_h2odata does not match value in GFS_typedefs.F90: ", & + h2o_coeff, " /= ", size(Data(1)%Tbd%h2opl, dim=3) + errflg = 1 + return + end if + + ! DH* OpenMP parallel region with OpenMP do loops for the three do loops? + ! there is no threading on the outside, i.e. can use CCPP's omp_threads here + + ! Update values of oz_pres in Interstitial data type for all threads + if (Model%ntoz > 0) then + do nt=1,nthrds + Interstitial(nt)%oz_pres = oz_pres + end do + end if + + ! Update values of h2o_pres in Interstitial data type for all threads + if (Model%h2o_phys) then + do nt=1,nthrds + Interstitial(nt)%h2o_pres = h2o_pres + end do + end if + + !--- read in and initialize ozone and water + if (Model%ntoz > 0) then + do nb = 1, nblks + call setindxoz (Model%blksz(nb), Data(nb)%Grid%xlat_d, Data(nb)%Grid%jindx1_o3, & + Data(nb)%Grid%jindx2_o3, Data(nb)%Grid%ddy_o3) + enddo + endif + + if (Model%h2o_phys) then + do nb = 1, nblks + call setindxh2o (Model%blksz(nb), Data(nb)%Grid%xlat_d, Data(nb)%Grid%jindx1_h, & + Data(nb)%Grid%jindx2_h, Data(nb)%Grid%ddy_h) + enddo + endif + + ! *DH + + is_initialized = .true. + end subroutine GFS_phys_time_vary_init - subroutine GFS_phys_time_vary_finalize() + +!> \section arg_table_GFS_phys_time_vary_finalize Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! + subroutine GFS_phys_time_vary_finalize(errmsg, errflg) + + implicit none + + ! Interface variables + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (.not.is_initialized) return + + ! Deallocate ozone arrays + if (allocated(oz_lat) ) deallocate(oz_lat) + if (allocated(oz_pres) ) deallocate(oz_pres) + if (allocated(oz_time) ) deallocate(oz_time) + if (allocated(ozplin) ) deallocate(ozplin) + + ! Deallocate h2o arrays + if (allocated(h2o_lat) ) deallocate(h2o_lat) + if (allocated(h2o_pres)) deallocate(h2o_pres) + if (allocated(h2o_time)) deallocate(h2o_time) + if (allocated(h2oplin) ) deallocate(h2oplin) + + is_initialized = .false. + end subroutine GFS_phys_time_vary_finalize + !> \section arg_table_GFS_phys_time_vary_run Argument Table !! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | !! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| @@ -25,9 +179,8 @@ end subroutine GFS_phys_time_vary_finalize !! subroutine GFS_phys_time_vary_run (Data, Model, errmsg, errflg) - use mersenne_twister, only: random_setseed, random_number + use mersenne_twister, only: random_setseed, random_number use machine, only: kind_phys - use physcons, only: dxmin, dxinv use GFS_typedefs, only: GFS_control_type, GFS_data_type implicit none @@ -52,6 +205,13 @@ subroutine GFS_phys_time_vary_run (Data, Model, errmsg, errflg) errmsg = '' errflg = 0 + ! Check initialization status + if (.not.is_initialized) then + write(errmsg,'(*(a))') "Logic error: GFS_phys_time_vary_run called before GFS_phys_time_vary_init" + errflg = 1 + return + end if + nblks = size(Model%blksz) !--- switch for saving convective clouds - cnvc90.f @@ -141,4 +301,4 @@ subroutine GFS_phys_time_vary_run (Data, Model, errmsg, errflg) end subroutine GFS_phys_time_vary_run - end module GFS_phys_time_vary + end module GFS_phys_time_vary diff --git a/physics/GFS_stochastics.F90 b/physics/GFS_stochastics.F90 index 390bffd5a..70075fc12 100644 --- a/physics/GFS_stochastics.F90 +++ b/physics/GFS_stochastics.F90 @@ -44,6 +44,8 @@ end subroutine GFS_stochastics_finalize !! | tprcp | nonnegative_lwe_thickness_of_precipitation_amount_on_dynamics_timestep | total precipitation amount in each time step | m | 1 | real | kind_phys | inout | F | !! | totprcp | accumulated_lwe_thickness_of_precipitation_amount | accumulated total precipitation | m | 1 | real | kind_phys | inout | F | !! | cnvprcp | cumulative_lwe_thickness_of_convective_precipitation_amount | cumulative convective precipitation | m | 1 | real | kind_phys | inout | F | +!! | totprcpb | accumulated_lwe_thickness_of_precipitation_amount_in_bucket | accumulated total precipitation in bucket | m | 1 | real | kind_phys | inout | F | +!! | cnvprcpb | cumulative_lwe_thickness_of_convective_precipitation_amount_in_bucket | cumulative convective precipitation in bucket | m | 1 | real | kind_phys | inout | F | !! | cplflx | flag_for_flux_coupling | flag controlling cplflx collection (default off) | flag | 0 | logical | | in | F | !! | rain_cpl | lwe_thickness_of_precipitation_amount_for_coupling | total rain precipitation | m | 1 | real | kind_phys | inout | F | !! | snow_cpl | lwe_thickness_of_snow_amount_for_coupling | total snow precipitation | m | 1 | real | kind_phys | inout | F | @@ -68,7 +70,8 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, shum_wts_inv, diss_est, & ugrs, vgrs, tgrs, qgrs, gu0, gv0, gt0, gq0, dtdtr, & rain, rainc, tprcp, totprcp, cnvprcp, & - cplflx, rain_cpl, snow_cpl, drain_cpl, dsnow_cpl, & + totprcpb, cnvprcpb, cplflx, & + rain_cpl, snow_cpl, drain_cpl, dsnow_cpl, & errmsg, errflg) use machine, only: kind_phys @@ -110,6 +113,8 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, real(kind_phys), dimension(1:im), intent(inout) :: tprcp real(kind_phys), dimension(1:im), intent(inout) :: totprcp real(kind_phys), dimension(1:im), intent(inout) :: cnvprcp + real(kind_phys), dimension(1:im), intent(inout) :: totprcpb + real(kind_phys), dimension(1:im), intent(inout) :: cnvprcpb logical, intent(in) :: cplflx ! rain_cpl, snow_cpl only allocated if cplflx == .true. or do_sppt == .true. real(kind_phys), dimension(:), intent(inout) :: rain_cpl @@ -174,6 +179,10 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, totprcp(:) = totprcp(:) + (sppt_wts(:,15) - 1 )*rain(:) ! acccumulated total and convective preciptiation cnvprcp(:) = cnvprcp(:) + (sppt_wts(:,15) - 1 )*rainc(:) + ! bucket precipitation adjustment due to sppt + totprcpb(:) = totprcpb(:) + (sppt_wts(:,15) - 1 )*rain(:) + cnvprcpb(:) = cnvprcpb(:) + (sppt_wts(:,15) - 1 )*rainc(:) + if (cplflx) then rain_cpl(:) = rain_cpl(:) + (sppt_wts(:,15) - 1.0)*drain_cpl(:) snow_cpl(:) = snow_cpl(:) + (sppt_wts(:,15) - 1.0)*dsnow_cpl(:) diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index a4d80dc07..8d058e46f 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -105,6 +105,8 @@ end subroutine GFS_suite_interstitial_1_finalize !! | frland | land_area_fraction | land area fraction | frac | 1 | real | kind_phys | out | F | !! | work1 | grid_size_related_coefficient_used_in_scale-sensitive_schemes | grid size related coefficient used in scale-sensitive schemes | none | 1 | real | kind_phys | out | F | !! | work2 | grid_size_related_coefficient_used_in_scale-sensitive_schemes_complement | complement to work1 | none | 1 | real | kind_phys | out | F | +!! | dxmin | minimum_scaling_factor_for_critical_relative_humidity | minimum scaling factor for critical relative humidity | m2 rad-2 | 0 | real | kind_phys | in | F | +!! | dxinv | inverse_scaling_factor_for_critical_relative_humidity | inverse scaling factor for critical relative humidity | rad2 m-2 | 0 | real | kind_phys | in | F | !! | dudt | tendency_of_x_wind_due_to_model_physics | updated tendency of the x wind | m s-2 | 2 | real | kind_phys | out | F | !! | dvdt | tendency_of_y_wind_due_to_model_physics | updated tendency of the y wind | m s-2 | 2 | real | kind_phys | out | F | !! | dtdt | tendency_of_air_temperature_due_to_model_physics | updated tendency of the temperature | K s-1 | 2 | real | kind_phys | out | F | @@ -114,10 +116,9 @@ end subroutine GFS_suite_interstitial_1_finalize !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! subroutine GFS_suite_interstitial_1_run (Model, Grid, Sfcprop, Statein, Diag, rhbbot, rhpbl, rhbtop, frain, islmsk, & - frland, work1, work2, dudt, dvdt, dtdt, dtdtc, dqdt, errmsg, errflg) + frland, work1, work2, dxmin, dxinv, dudt, dvdt, dtdt, dtdtc, dqdt, errmsg, errflg) use machine, only: kind_phys - use physcons, only: dxmin, dxinv use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_sfcprop_type, GFS_statein_type, GFS_diag_type implicit none @@ -135,6 +136,7 @@ subroutine GFS_suite_interstitial_1_run (Model, Grid, Sfcprop, Statein, Diag, rh real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(out) :: work1, work2 real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(out) :: dudt, dvdt, dtdt, dtdtc real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs,Model%ntrac), intent(out) :: dqdt + real(kind=kind_phys), intent(in) :: dxmin, dxinv character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -370,15 +372,15 @@ end subroutine GFS_suite_interstitial_3_finalize !! | ktop | vertical_index_at_cloud_top | vertical index at cloud top | index | 1 | integer | | inout | F | !! | kbot | vertical_index_at_cloud_base | vertical index at cloud base | index | 1 | integer | | inout | F | !! | rhc | critical_relative_humidity | critical relative humidity | frac | 2 | real | kind_phys | out | F | +!! | rhcmax | maximum_critical_relative_humidity | maximum critical relative humidity | frac | 0 | real | kind_phys | in | F | !! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! subroutine GFS_suite_interstitial_3_run (Model, Grid, Statein, rhbbot, rhbtop, work1, work2, clw, cnvc, cnvw, & - ktop, kbot, rhc, errmsg, errflg) + ktop, kbot, rhc, rhcmax, errmsg, errflg) use machine, only: kind_phys use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_statein_type - use physcons, only: rhc_max implicit none @@ -392,6 +394,7 @@ subroutine GFS_suite_interstitial_3_run (Model, Grid, Statein, rhbbot, rhbtop, w real(kind=kind_phys), intent(inout) :: clw(:,:,:), cnvc(:,:), cnvw(:,:) integer, dimension(size(Grid%xlon,1)), intent(out) :: ktop, kbot real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(out) :: rhc + real(kind=kind_phys), intent(in) :: rhcmax character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -419,7 +422,7 @@ subroutine GFS_suite_interstitial_3_run (Model, Grid, Statein, rhbbot, rhbtop, w do k=1,Model%levs do i=1, size(Grid%xlon,1) tem = rhbbot - (rhbbot-rhbtop) * (1.0-Statein%prslk(i,k)) - tem = rhc_max * work1(i) + tem * work2(i) + tem = rhcmax * work1(i) + tem * work2(i) rhc(i,k) = max(0.0, min(1.0,tem)) enddo enddo diff --git a/physics/GFS_time_vary_pre.fv3.f90 b/physics/GFS_time_vary_pre.fv3.f90 index 11b8b7119..d0c604879 100644 --- a/physics/GFS_time_vary_pre.fv3.f90 +++ b/physics/GFS_time_vary_pre.fv3.f90 @@ -3,20 +3,63 @@ module GFS_time_vary_pre + use funcphys, only: gfuncphys + implicit none private public GFS_time_vary_pre_init, GFS_time_vary_pre_run, GFS_time_vary_pre_finalize + logical :: is_initialized = .false. + contains - subroutine GFS_time_vary_pre_init () +!> \section arg_table_GFS_time_vary_pre_init Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! + subroutine GFS_time_vary_pre_init (errmsg, errflg) + + implicit none + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + if (is_initialized) return + + !--- Call gfuncphys (funcphys.f) to compute all physics function tables. + call gfuncphys () + + is_initialized = .true. + end subroutine GFS_time_vary_pre_init - subroutine GFS_time_vary_pre_finalize() + +!> \section arg_table_GFS_time_vary_pre_finalize Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! + subroutine GFS_time_vary_pre_finalize(errmsg, errflg) + + implicit none + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + if (.not. is_initialized) return + + ! DH* this is the place to deallocate whatever is allocated by gfuncphys() in GFS_time_vary_pre_init + + is_initialized = .false. + end subroutine GFS_time_vary_pre_finalize + !> \section arg_table_GFS_time_vary_pre_run Argument Table !! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | !! |----------------|--------------------------------------------------------|-------------------------------------------------------------------------|----------|------|-----------------------|-----------|--------|----------| @@ -43,6 +86,13 @@ subroutine GFS_time_vary_pre_run (Model, errmsg, errflg) errmsg = '' errflg = 0 + ! Check initialization status + if (.not.is_initialized) then + write(errmsg,'(*(a))') "Logic error: GFS_time_vary_pre_run called before GFS_time_vary_pre_init" + errflg = 1 + return + end if + !--- Model%jdat is being updated directly inside of FV3GFS_cap.F90 !--- update calendars and triggers rinc(1:5) = 0 diff --git a/physics/gfdl_cloud_microphys.F90 b/physics/gfdl_cloud_microphys.F90 index d56975ed4..9061604c9 100644 --- a/physics/gfdl_cloud_microphys.F90 +++ b/physics/gfdl_cloud_microphys.F90 @@ -51,9 +51,12 @@ module gfdl_cloud_microphys real :: missing_value = - 1.e10 character (len = 17) :: mod_name = 'gfdl_cloud_microphys' - ! DH* CLEANUP!!! + ! DH* TODO: CLEANUP, all of these should be coming in through the argument list real(kind=kind_phys), parameter :: one = 1.0d0 real(kind=kind_phys), parameter :: con_d00 = 0.0d0 + real(kind=kind_phys), parameter :: con_p001= 0.001d0 + real(kind=kind_phys), parameter :: con_day = 86400.d0 + real(kind=kind_phys), parameter :: rainmin = 1.0e-13 real, parameter :: grav = 9.80665 !< gfs: acceleration due to gravity real, parameter :: rgrav = 1./grav real, parameter :: rdgas = 287.05 !< gfs: gas constant for dry air @@ -368,7 +371,7 @@ subroutine gfdl_cloud_microphys_init (me, master, nlunit, input_nml_file, loguni if (is_initialized) return if (imp_physics/=imp_physics_gfdl) then - write(errmsg,'(*(a))') 'Namelist option for micrphysics does not match choice in suite definition file' + write(errmsg,'(*(a))') 'Namelist option for microphysics does not match choice in suite definition file' errflg = 1 return end if @@ -525,6 +528,7 @@ subroutine gfdl_cloud_microphys_run( & qs1, pt_dt, qa_dt, u_dt, v_dt, w, qv_dt, ql_dt, qr_dt, qi_dt, & qs_dt, qg_dt real(kind=kind_phys) :: onebg + real(kind=kind_phys) :: tem ! Initialize CCPP error handling variables errmsg = '' @@ -587,12 +591,26 @@ subroutine gfdl_cloud_microphys_run( & garea, dtp, frland, rain0, snow0, ice0, graupel0, hydrostatic, & phys_hydrostatic) + tem = dtp*con_p001/con_day + ! fix negative values do i = 1, im - rain0(i) = max(con_d00, rain0(i)) - snow0(i) = max(con_d00, snow0(i)) - ice0(i) = max(con_d00, ice0(i)) - graupel0(i) = max(con_d00, graupel0(i)) + !rain0(i) = max(con_d00, rain0(i)) + !snow0(i) = max(con_d00, snow0(i)) + !ice0(i) = max(con_d00, ice0(i)) + !graupel0(i) = max(con_d00, graupel0(i)) + if(rain0(i)*tem < rainmin) then + rain0(i) = 0.0 + endif + if(ice0(i)*tem < rainmin) then + ice0(i) = 0.0 + endif + if(snow0(i)*tem < rainmin) then + snow0(i) = 0.0 + endif + if(graupel0(i)*tem < rainmin) then + graupel0(i) = 0.0 + endif enddo ! flip vertical coordinate back diff --git a/physics/gfdl_cloud_microphys_pre_post.F90 b/physics/gfdl_cloud_microphys_pre_post.F90 index 655d33280..2a5a1f4b0 100644 --- a/physics/gfdl_cloud_microphys_pre_post.F90 +++ b/physics/gfdl_cloud_microphys_pre_post.F90 @@ -63,7 +63,7 @@ module gfdl_cloud_microphys_post public gfdl_cloud_microphys_post_run, gfdl_cloud_microphys_post_init, gfdl_cloud_microphys_post_finalize - ! DH* CLEANUP !!! + ! DH* TODO: CLEANUP, all of these should be coming in through the argument list real(kind=kind_phys), parameter :: con_p001= 0.001d0 real(kind=kind_phys), parameter :: con_day = 86400.d0 real(kind=kind_phys), parameter :: rainmin = 1.0e-13 diff --git a/physics/h2ointerp.f90 b/physics/h2ointerp.f90 index af3021f4e..82a765e71 100644 --- a/physics/h2ointerp.f90 +++ b/physics/h2ointerp.f90 @@ -1,3 +1,13 @@ +module h2ointerp + + implicit none + + private + + public :: read_h2odata, setindxh2o, h2ointerpol + +contains + subroutine read_h2odata (h2o_phys, me, master) use machine, only: kind_phys use h2o_def @@ -185,3 +195,5 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) ! return end + +end module h2ointerp \ No newline at end of file diff --git a/physics/h2ophys.f b/physics/h2ophys.f new file mode 100755 index 000000000..8ee3c375c --- /dev/null +++ b/physics/h2ophys.f @@ -0,0 +1,157 @@ + module h2ophys + + implicit none + + private + + public :: h2ophys_init, h2ophys_run, h2ophys_finalize + + contains + +! \brief Brief description of the subroutine +! +!> \section arg_table_h2ophys_init Argument Table +!! + subroutine h2ophys_init() + end subroutine h2ophys_init + +!>\defgroup GFS_h2ophys GFS h2ophys Main +!! \section arg_table_h2ophys_run Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|---------------------------------------------------|---------------------------------------------------|---------|------|-----------|-----------|--------|----------| +!! | ix | horizontal_dimension | horizontal dimension | count | 0 | integer | | in | F | +!! | im | horizontal_loop_extent | horizontal loop extent | count | 0 | integer | | in | F | +!! | levs | vertical_dimension | number of vertical layers | count | 0 | integer | | in | F | +!! | kh2o | vertical_dimension_of_h2o_forcing_data | number of vertical layers in h2o forcing data | count | 0 | integer | | in | F | +!! | dt | time_step_for_physics | physics time step | s | 0 | real | kind_phys | in | F | +!! | h2o | water_vapor_specific_humidity_updated_by_physics | water vapor specific humidity updated by physics | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | ph2o | natural_log_of_h2o_forcing_data_pressure_levels | natural log of h2o forcing data pressure levels | log(Pa) | 1 | real | kind_phys | in | F | +!! | prsl | air_pressure | mid-layer pressure | Pa | 2 | real | kind_phys | in | F | +!! | h2opltc | h2o_forcing | water forcing data | various | 3 | real | kind_phys | in | F | +!! | h2o_coeff | number_of_coefficients_in_h2o_forcing_data | number of coefficients in h2o forcing data | index | 0 | integer | | in | F | +!! | ldiag3d | flag_diagnostics_3D | flag for calculating 3-D diagnostic fields | flag | 0 | logical | | in | F | +!! | h2op | change_in_h2o_concentration | change in h2o concentration | kg kg-1 | 3 | real | kind_phys | inout | F | +!! | me | mpi_rank | rank of the current MPI task | index | 0 | integer | | in | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +!! \section genal_h2ophys GFS H2O Physics Scheme General Algorithm +!> @{ + subroutine h2ophys_run(ix, im, levs, kh2o, dt, h2o, ph2o, prsl, + & h2opltc, h2o_coeff, ldiag3d, h2op, me, + & errmsg, errflg) +! +! May 2015 - Shrinivas Moorthi - Adaptation of NRL H2O physics for +! stratosphere and mesosphere +! +! this code assumes that both prsl and ph2o are from bottom to top +! as are all other variables +! + use machine , only : kind_phys + implicit none +! interface variables + integer, intent(in) :: ix, im, levs, kh2o, h2o_coeff, me + real(kind=kind_phys), intent(in) :: dt + real(kind=kind_phys), intent(inout) :: h2o(ix,levs) + real(kind=kind_phys), intent(in) :: ph2o(kh2o) + real(kind=kind_phys), intent(in) :: prsl(ix,levs) + real(kind=kind_phys), intent(in) :: h2opltc(ix,kh2o,h2o_coeff) + logical , intent(in) :: ldiag3d + real(kind=kind_phys), intent(inout) :: h2op(ix,levs,h2o_coeff) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg +! local variables + integer k,kmax,kmin,l,i,j + logical flg(im) + real(kind=kind_phys) pmax, pmin, tem, temp + real(kind=kind_phys) wk1(im), wk2(im), wk3(im), pltc(im,h2o_coeff) + &, h2oib(im) + real, parameter :: prsmax=10000.0, pmaxl=log(prsmax) +! +! initialize CCPP error handling variables + errmsg = '' + errflg = 0 +! +! write(1000+me,*)' in h2ophys ix=',ix, im, levs, kh2o, dt + do l=1,levs + pmin = 1.0e10 + pmax = -1.0e10 +! + do i=1,im + wk1(i) = log(prsl(i,l)) + pmin = min(wk1(i), pmin) + pmax = max(wk1(i), pmax) + pltc(i,:) = 0.0 + enddo + if (pmin < pmaxl) then + kmax = 1 + kmin = 1 + do k=1,kh2o-1 + if (pmin < ph2o(k)) kmax = k + if (pmax < ph2o(k)) kmin = k + enddo +! + do k=kmin,kmax + temp = 1.0 / (ph2o(k) - ph2o(k+1)) + do i=1,im + flg(i) = .false. + if (wk1(i) < ph2o(k) .and. wk1(i) >= ph2o(k+1)) then + flg(i) = .true. + wk2(i) = (wk1(i) - ph2o(k+1)) * temp + wk3(i) = 1.0 - wk2(i) + endif + enddo + do j=1,h2o_coeff + do i=1,im + if (flg(i)) then + pltc(i,j) = wk2(i) * h2opltc(i,k,j) + & + wk3(i) * h2opltc(i,k+1,j) + endif + enddo + enddo + enddo +! + do j=1,h2o_coeff + do i=1,im + if (wk1(i) < ph2o(kh2o)) then + pltc(i,j) = h2opltc(i,kh2o,j) + endif + if (wk1(i) >= ph2o(1)) then + pltc(i,j) = h2opltc(i,1,j) + endif + enddo + enddo + endif + do i=1,im + if (prsl(i,l) < prsmax) then + h2oib(i) = h2o(i,l) ! no filling + tem = 1.0 / pltc(i,2) ! 1/teff + h2o(i,l) = (h2oib(i) + (pltc(i,1)+pltc(i,3)*tem)*dt) + & / (1.0 + tem*dt) + endif + +! if (i == 1) write(1000+me,*)' h2oib=',h2oib(i),' pltc1=', +! &pltc(i,1),' pltc2=', pltc(i,2),' tem=',tem ,' dt=',dt +! &,' l=',l + enddo +! +! if (ldiag3d) then ! h2o change diagnostics +! do i=1,im +! h2op(i,l,1) = h2op(i,l,1) + pltc(i,1)*dt +! h2op(i,l,2) = h2op(i,l,2) + (h2oo(i,l) - h2oib(i)) +! enddo +! endif + enddo ! vertical loop +! + return + end subroutine h2ophys_run +!> @} + +! \brief Brief description of the subroutine +! +!> \section arg_table_h2ophys_finalize Argument Table +!! + subroutine h2ophys_finalize() + end subroutine h2ophys_finalize + + end module h2ophys \ No newline at end of file diff --git a/physics/ozinterp.f90 b/physics/ozinterp.f90 index 656bfafbe..e7704924e 100644 --- a/physics/ozinterp.f90 +++ b/physics/ozinterp.f90 @@ -1,3 +1,13 @@ +module ozinterp + + implicit none + + private + + public :: read_o3data, setindxoz, ozinterpol + +contains + SUBROUTINE read_o3data (ntoz, me, master) use machine, only: kind_phys use ozne_def @@ -9,6 +19,7 @@ SUBROUTINE read_o3data (ntoz, me, master) integer :: i, n, k real(kind=4), allocatable, dimension(:) :: oz_lat4, oz_pres4 real(kind=4), allocatable, dimension(:) :: oz_time4, tempin + real(kind=4) :: blatc4 if (ntoz <= 0) then ! Diagnostic ozone rewind (kozc) @@ -191,3 +202,5 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) ! RETURN END + +end module ozinterp \ No newline at end of file diff --git a/physics/ozphys.f b/physics/ozphys.f index b872a0dcc..f5e276589 100644 --- a/physics/ozphys.f +++ b/physics/ozphys.f @@ -23,7 +23,6 @@ end subroutine ozphys_pre_run subroutine ozphys_pre_finalize() end subroutine ozphys_pre_finalize - end module ozphys_pre @@ -56,7 +55,7 @@ end subroutine ozphys_init !! | po3 | natural_log_of_ozone_forcing_data_pressure_levels | natural log of ozone forcing data pressure levels | log(Pa) | 1 | real | kind_phys | in | F | !! | prsl | air_pressure | mid-layer pressure | Pa | 2 | real | kind_phys | in | F | !! | prdout | ozone_forcing | ozone forcing data | various | 3 | real | kind_phys | in | F | -!! | pl_coeff | number_of_coefficients_in_ozone_forcing_data | number of coefficients in ozone forcing data | index | 0 | integer | | in | F | +!! | oz_coeff | number_of_coefficients_in_ozone_forcing_data | number of coefficients in ozone forcing data | index | 0 | integer | | in | F | !! | delp | air_pressure_difference_between_midlayers | difference between mid-layer pressures | Pa | 2 | real | kind_phys | in | F | !! | ldiag3d | flag_diagnostics_3D | flag for calculating 3-D diagnostic fields | flag | 0 | logical | | in | F | !! | ozp | change_in_ozone_concentration | change in ozone concentration | kg kg-1 | 3 | real | kind_phys | inout | F | @@ -68,7 +67,7 @@ end subroutine ozphys_init !> @{ subroutine ozphys_run ( & & ix, im, levs, ko3, dt, oz, tin, po3, & - & prsl, prdout, pl_coeff, delp, ldiag3d, & + & prsl, prdout, oz_coeff, delp, ldiag3d, & & ozp, me, errmsg, errflg) ! ! this code assumes that both prsl and po3 are from bottom to top @@ -79,11 +78,11 @@ subroutine ozphys_run ( & implicit none ! ! Interface variables - integer, intent(in) :: im, ix, levs, ko3, pl_coeff, me + integer, intent(in) :: im, ix, levs, ko3, oz_coeff, me real(kind=kind_phys), intent(inout) :: & - & oz(ix,levs), ozp(ix,levs,pl_coeff) + & oz(ix,levs), ozp(ix,levs,oz_coeff) real(kind=kind_phys), intent(in) :: & - & dt, po3(ko3), prdout(ix,ko3,pl_coeff), & + & dt, po3(ko3), prdout(ix,ko3,oz_coeff), & & prsl(ix,levs), tin(ix,levs), delp(ix,levs) logical, intent(in) :: ldiag3d character(len=*), intent(out) :: errmsg @@ -94,7 +93,7 @@ subroutine ozphys_run ( & integer k,kmax,kmin,l,i,j logical flg(im) real(kind=kind_phys) pmax, pmin, tem, temp - real(kind=kind_phys) wk1(im), wk2(im), wk3(im), prod(im,pl_coeff), + real(kind=kind_phys) wk1(im), wk2(im), wk3(im), prod(im,oz_coeff), & ozib(im), colo3(im,levs+1), ozi(ix,levs) ! ! Initialize CCPP error handling variables @@ -105,7 +104,7 @@ subroutine ozphys_run ( & ozi = oz ! !> - Calculate vertical integrated column ozone values. - if (pl_coeff > 2) then + if (oz_coeff > 2) then colo3(:,levs+1) = 0.0 do l=levs,1,-1 do i=1,im @@ -142,7 +141,7 @@ subroutine ozphys_run ( & wk3(i) = 1.0 - wk2(i) endif enddo - do j=1,pl_coeff + do j=1,oz_coeff do i=1,im if (flg(i)) then prod(i,j) = wk2(i) * prdout(i,k,j) @@ -152,7 +151,7 @@ subroutine ozphys_run ( & enddo enddo ! - do j=1,pl_coeff + do j=1,oz_coeff do i=1,im if (wk1(i) < po3(ko3)) then prod(i,j) = prdout(i,ko3,j) @@ -163,7 +162,7 @@ subroutine ozphys_run ( & enddo enddo - if (pl_coeff == 2) then + if (oz_coeff == 2) then do i=1,im ozib(i) = ozi(i,l) ! no filling oz(i,l) = (ozib(i) + prod(i,1)*dt) / (1.0 + prod(i,2)*dt) @@ -181,7 +180,7 @@ subroutine ozphys_run ( & !! - ozp(:,:,2) - Ozone tendency at model layers !! - ozp(:,:,3) - Ozone production from temperature term at model layers !! - ozp(:,:,4) - Ozone production from column ozone term at model layers - if (pl_coeff == 4) then + if (oz_coeff == 4) then do i=1,im ozib(i) = ozi(i,l) ! no filling tem = prod(i,1) + prod(i,3)*tin(i,l) @@ -226,31 +225,34 @@ subroutine ozphys_post_init() end subroutine ozphys_post_init +#if 0 +! DH* TODO - split dq3dt into individual components? Or replace Interstitial%dq3dt_loc(:,:,1:) / dq3dt_loc(:,:,1:) +! in GFS_physics_driver entirely with Diag%dq3dt(:,:,6:) ? *DH !! \section arg_table_ozphys_post_run Argument Table -!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -!! |----------------|----------------------------------------------|----------------------------------------------|---------|------|---------------|-----------|--------|----------| -!! | ix | horizontal_dimension | horizontal dimension | count | 0 | integer | | in | F | -!! | levs | vertical_dimension | number of vertical layers | count | 0 | integer | | in | F | -!! | pl_coeff | number_of_coefficients_in_ozone_forcing_data | number of coefficients in ozone forcing data | index | 0 | integer | | in | F | -!! | ldiag3d | flag_diagnostics_3D | logical flag for 3D diagnostics | flag | 0 | logical | | in | F | -!! | ozp | change_in_ozone_concentration | change in ozone concentration | kg kg-1 | 3 | real | kind_phys | in | F | -!! | Diag | FV3-GFS_Diag_type | GFS diagnostics derived data type variable | DDT | 0 | GFS_diag_type | | inout | F | -!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | -!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|-------------------------------------------------------------------|-------------------------------------------------------------------|---------|------|---------------|-----------|--------|----------| +!! | im | horizontal_loop_extent | horizontal loop extent | count | 0 | integer | | in | F | +!! | levs | vertical_dimension | number of vertical layers | count | 0 | integer | | in | F | +!! | oz_coeff | number_of_coefficients_in_ozone_forcing_data | number of coefficients in ozone forcing data | index | 0 | integer | | in | F | +!! | ldiag3d | flag_diagnostics_3D | logical flag for 3D diagnostics | flag | 0 | logical | | in | F | +!! | ozp | change_in_ozone_concentration | change in ozone concentration | kg kg-1 | 3 | real | kind_phys | in | F | +!! | dq3dt | cumulative_change_in_water_vapor_specific_humidity_due_to_physics | cumulative change in water vapor specific humidity due to physics | kg kg-1 | 3 | real | kind_phys | inout | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! - subroutine ozphys_post_run(ix, levs, pl_coeff, ldiag3d, ozp, & - & Diag, errmsg, errflg) +#endif + subroutine ozphys_post_run(im, levs, oz_coeff, ldiag3d, ozp, & + & dq3dt, errmsg, errflg) use GFS_typedefs, only: GFS_diag_type use machine, only: kind_phys implicit none - integer, intent(in) :: ix, levs, pl_coeff + integer, intent(in) :: im, levs, oz_coeff logical, intent(in) :: ldiag3d - real(kind=kind_phys), intent(in) :: ozp(ix,levs,pl_coeff) - type(GFS_diag_type), intent(inout) :: Diag - + real(kind=kind_phys), intent(in) :: ozp(im,levs,oz_coeff) + real(kind=kind_phys), intent(inout) :: dq3dt(im,levs,oz_coeff+5) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -259,10 +261,10 @@ subroutine ozphys_post_run(ix, levs, pl_coeff, ldiag3d, ozp, & errflg = 0 if (ldiag3d) then - Diag%dq3dt(:,:,6) = ozp(:,:,1) - Diag%dq3dt(:,:,7) = ozp(:,:,2) - Diag%dq3dt(:,:,8) = ozp(:,:,3) - Diag%dq3dt(:,:,9) = ozp(:,:,4) + dq3dt(:,:,6) = ozp(:,:,1) + dq3dt(:,:,7) = ozp(:,:,2) + dq3dt(:,:,8) = ozp(:,:,3) + dq3dt(:,:,9) = ozp(:,:,4) end if return diff --git a/physics/physcons.f90 b/physics/physcons.F90 similarity index 99% rename from physics/physcons.f90 rename to physics/physcons.F90 index 6e1d6ae00..e94e604e2 100644 --- a/physics/physcons.f90 +++ b/physics/physcons.F90 @@ -178,8 +178,6 @@ module physcons ! real(kind=kind_phys),parameter:: rhosnow = 100. ! density of snow (kg/m^3) real(kind=kind_phys),parameter:: rhoair = 1.28 ! density of air near surface (kg/m^3) - real(kind=kind_phys) :: dxmax, dxmin, dxinv, rhc_max - !........................................! end module physcons ! !========================================! diff --git a/physics/set_soilveg.f b/physics/set_soilveg.f index 333748b9d..17f9e4ae9 100644 --- a/physics/set_soilveg.f +++ b/physics/set_soilveg.f @@ -1,9 +1,18 @@ + module set_soilveg_mod + + implicit none + + private + + public set_soilveg + + contains + subroutine set_soilveg(me,isot,ivet,nlunit) use namelist_soilveg implicit none - integer, intent(in) :: isot,ivet,nlunit - integer me + integer, intent(in) :: me,isot,ivet,nlunit !my begin locals !for 20 igbp veg type and 19 stasgo soil type integer i @@ -405,4 +414,6 @@ subroutine set_soilveg(me,isot,ivet,nlunit) ! if (me == 0) write(6,soil_veg) return - end + end subroutine set_soilveg + + end module set_soilveg_mod diff --git a/physics/sfc_drv.f b/physics/sfc_drv.f index 513aae26e..a2fb070dc 100644 --- a/physics/sfc_drv.f +++ b/physics/sfc_drv.f @@ -134,14 +134,76 @@ end module lsm_noah_post !! @} module lsm_noah + + use set_soilveg_mod, only: set_soilveg + + implicit none + + private + + public :: lsm_noah_init, lsm_noah_run, lsm_noah_finalize + + logical :: is_initialized = .false. + contains - subroutine lsm_noah_init +!! \section arg_table_lsm_noah_init Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|-------------------------------------------------------------|---------------------------------------------------------|------------|------|-----------|-----------|--------|----------| +!! | me | mpi_rank | current MPI-rank | index | 0 | integer | | in | F | +!! | isot | soil_type_dataset_choice | soil type dataset choice | index | 0 | integer | | in | F | +!! | ivegsrc | vegetation_type_dataset_choice | land use dataset choice | index | 0 | integer | | in | F | +!! | nlunit | iounit_namelist | fortran unit number for file opens | none | 0 | integer | | in | F | +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! + subroutine lsm_noah_init(me, isot, ivegsrc, nlunit, + & errmsg, errflg) + + implicit none + + integer, intent(in) :: me, isot, ivegsrc, nlunit + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (is_initialized) return + + !--- initialize soil vegetation + call set_soilveg(me, isot, ivegsrc, nlunit) + + is_initialized = .true. + end subroutine lsm_noah_init - subroutine lsm_noah_finalize + +!! \section arg_table_lsm_noah_finalize Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|-------------------------------------------------------------|--------------------------------------------|------------|------|-----------|-----------|--------|----------| +!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! + subroutine lsm_noah_finalize(errmsg, errflg) + + implicit none + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (.not. is_initialized) return + + is_initialized = .false. + end subroutine lsm_noah_finalize + ! ===================================================================== ! ! description: ! ! ! @@ -439,6 +501,14 @@ subroutine lsm_noah_run & errmsg = '' errflg = 0 + ! Check initialization status + if (.not. is_initialized) then + write(errmsg,'(*(a))') + & "Logic error: lsm_noah_run called before lsm_noah_init" + errflg = 1 + return + end if + !> - Set flag for land points. do i = 1, im diff --git a/physics/sflx.f b/physics/sflx.f index d9ae482ca..3f4105bfe 100644 --- a/physics/sflx.f +++ b/physics/sflx.f @@ -6,14 +6,14 @@ !! It is a soil/veg/snowpack land-surface model to update soil moisture, soil !! ice, soil temperature, skin temperature, snowpack water content, snowdepth, !! and all terms of the surface energy balance and surface water balance -!! (excluding input atmospheric forcings of downward radiation and +!! (excluding input atmospheric forcings of downward radiation and !! precipitation ). !! !! The land-surface model component was substantially upgraded from the Oregon !! State University (OSU) land surface model to EMC's new Noah Land Surface Model !! (Noah LSM) during the major implementation in the NCEP Global Forecast System !! (GFS) on May 31, 2005. Forecast System (GFS). The Noah LSM embodies about 10 -!! years of upgrades (see \cite chen_et_al_1996; \cite koren_et_al_1999; +!! years of upgrades (see \cite chen_et_al_1996; \cite koren_et_al_1999; !! \cite ek_et_al_2003) to its ancestor, the OSU LSM. The Noah LSM upgrade includes: !! - An increase from two (10, 190 cm thick) to four soil layers (10, 30, 60, 100 cm thick) !! - Addition of frozen soil physics @@ -28,83 +28,83 @@ !! - Improved seasonality of green vegetation cover. !! - Improved evaporation treatment over bare soil and snowpack !! -!!\param[in] nsoil integer, number of soil layers (>=2 but <=nsold) -!!\param[in] couple integer, =0:uncoupled (land model only), -!! =1:coupled with parent atmos model -!!\param[in] icein integer, sea-ice flag (=1: sea-ice, =0: land) -!!\param[in] ffrozp real, flag for snow-rain detection (1.=snow, 0.=rain) -!!\param[in] dt real, time step (<3600 sec) +!!\param[in] nsoil integer, number of soil layers (>=2 but <=nsold) +!!\param[in] couple integer, =0:uncoupled (land model only), +!! =1:coupled with parent atmos model +!!\param[in] icein integer, sea-ice flag (=1: sea-ice, =0: land) +!!\param[in] ffrozp real, flag for snow-rain detection (1.=snow, 0.=rain) +!!\param[in] dt real, time step (<3600 sec) !!\param[in] zlvl real, height abv atmos ground forcing vars (\f$m\f$) -!!\param[in] sldpth real, thickness of each soil layer (\f$m\f$), nsoil +!!\param[in] sldpth real, thickness of each soil layer (\f$m\f$), nsoil !!\param[in] swdn real, downward SW radiation flux (\f$W/m^2\f$) !!\param[in] swnet real, downward SW net (dn-up) flux (\f$W/m^2\f$) !!\param[in] lwdn real, downward LW radiation flux (\f$W/m^2\f$) -!!\param[in] sfcems real, sfc LW emissivity (fractional) -!!\param[in] sfcprs real, pressure at height zlvl above ground(\f$Pa\f$) -!!\param[in] sfctmp real, air temp at height zlvl above ground (\f$K\f$) -!!\param[in] sfcspd real, wind speed at height zlvl above ground (\f$m s^{-1}\f$) -!!\param[in] prcp real, precipitation rate (\f$kgm^{-2}s^{-1}\f$) -!!\param[in] q2 real, mixing ratio at hght zlvl above ground (\f$kgkg^{-1}\f$) -!!\param[in] q2sat real, sat mixing ratio at zlvl above ground (\f$kgkg^{-1}\f$) -!!\param[in] dqsdt2 real, slope of sat specific humidity curve at t=sfctmp (\f$kgkg^{-1}k^{-1}\f$) -!!\param[in] th2 real, air potential temperature at zlvl above ground (\f$K\f$) -!!\param[in] ivegsrc integer, sfc veg type data source UMD or IGBP -!!\param[in] vegtyp integer, vegetation type (integer index) -!!\param[in] soiltyp integer, soil type (integer index) -!!\param[in] slopetyp integer, class of sfc slope (integer index) -!!\param[in] shdmin real, min areal coverage of green veg (fraction) -!!\param[in] alb real, background snow-free sfc albedo (fraction) -!!\param[in] snoalb real, max albedo over deep snow (fraction) -!!\param[in] bexpp real, perturbation of soil type "b" parameter (perturbation) -!!\param[in] xlaip real, perturbation of leave area index (perturbation) -!!\param[in,out] tbot real, bottom soil temp (\f$K\f$) (local yearly-mean sfc air temp) +!!\param[in] sfcems real, sfc LW emissivity (fractional) +!!\param[in] sfcprs real, pressure at height zlvl above ground(\f$Pa\f$) +!!\param[in] sfctmp real, air temp at height zlvl above ground (\f$K\f$) +!!\param[in] sfcspd real, wind speed at height zlvl above ground (\f$m s^{-1}\f$) +!!\param[in] prcp real, precipitation rate (\f$kgm^{-2}s^{-1}\f$) +!!\param[in] q2 real, mixing ratio at hght zlvl above ground (\f$kgkg^{-1}\f$) +!!\param[in] q2sat real, sat mixing ratio at zlvl above ground (\f$kgkg^{-1}\f$) +!!\param[in] dqsdt2 real, slope of sat specific humidity curve at t=sfctmp (\f$kgkg^{-1}k^{-1}\f$) +!!\param[in] th2 real, air potential temperature at zlvl above ground (\f$K\f$) +!!\param[in] ivegsrc integer, sfc veg type data source UMD or IGBP +!!\param[in] vegtyp integer, vegetation type (integer index) +!!\param[in] soiltyp integer, soil type (integer index) +!!\param[in] slopetyp integer, class of sfc slope (integer index) +!!\param[in] shdmin real, min areal coverage of green veg (fraction) +!!\param[in] alb real, background snow-free sfc albedo (fraction) +!!\param[in] snoalb real, max albedo over deep snow (fraction) +!!\param[in] bexpp real, perturbation of soil type "b" parameter (perturbation) +!!\param[in] xlaip real, perturbation of leave area index (perturbation) +!!\param[in,out] tbot real, bottom soil temp (\f$K\f$) (local yearly-mean sfc air temp) !!\param[in,out] cmc real, canopy moisture content (\f$m\f$) -!!\param[in,out] t1 real, ground/canopy/snowpack eff skin temp (\f$K\f$) -!!\param[in,out] stc real, soil temp (\f$K\f$) -!!\param[in,out] smc real, total soil moisture (vol fraction) -!!\param[in,out] sh2o real, unfrozen soil moisture (vol fraction), note: frozen part = smc-sh2o -!!\param[in,out] sneqv real, water-equivalent snow depth (\f$m\f$), note: snow density = snwqv/snowh -!!\param[in,out] ch real, sfc exchange coeff for heat & moisture (\f$ms^{-1}\f$), note: conductance since it's been mult by wind +!!\param[in,out] t1 real, ground/canopy/snowpack eff skin temp (\f$K\f$) +!!\param[in,out] stc real, soil temp (\f$K\f$) +!!\param[in,out] smc real, total soil moisture (vol fraction) +!!\param[in,out] sh2o real, unfrozen soil moisture (vol fraction), note: frozen part = smc-sh2o +!!\param[in,out] sneqv real, water-equivalent snow depth (\f$m\f$), note: snow density = snwqv/snowh +!!\param[in,out] ch real, sfc exchange coeff for heat & moisture (\f$ms^{-1}\f$), note: conductance since it's been mult by wind !!\param[in,out] cm real, sfc exchange coeff for momentum (\f$ms^{-1}\f$), note: conductance since it's been mult by wind !!\param[in,out] z0 real, roughness length (\f$m\f$) -!!\param[out] nroot integer, number of root layers -!!\param[out] shdfac real, aeral coverage of green veg (fraction) -!!\param[out] snowh real, snow depth (\f$m\f$) -!!\param[out] albedo real, sfc albedo incl snow effect (fraction) -!!\param[out] eta real, downward latent heat flux (\f$W/m^2\f$) -!!\param[out] sheat real, downward sensible heat flux (\f$W/m^2\f$) -!!\param[out] ec real, canopy water evaporation (\f$W/m^2\f$) +!!\param[out] nroot integer, number of root layers +!!\param[out] shdfac real, aeral coverage of green veg (fraction) +!!\param[out] snowh real, snow depth (\f$m\f$) +!!\param[out] albedo real, sfc albedo incl snow effect (fraction) +!!\param[out] eta real, downward latent heat flux (\f$W/m^2\f$) +!!\param[out] sheat real, downward sensible heat flux (\f$W/m^2\f$) +!!\param[out] ec real, canopy water evaporation (\f$W/m^2\f$) !!\param[out] edir real, direct soil evaporation (\f$W/m^2\f$) !!\param[out] et real, plant transpiration (\f$W/m^2\f$) -!!\param[out] ett real, total plant transpiration (\f$W/m^2\f$) +!!\param[out] ett real, total plant transpiration (\f$W/m^2\f$) !!\param[out] esnow real, sublimation from snowpack (\f$W/m^2\f$) -!!\param[out] drip real, through-fall of precip and/or dew in excess of canopy water-holding capacity (\f$m\f$) -!!\param[out] dew real, dewfall (or frostfall for t<273.15) (\f$m\f$) -!!\param[out] beta real, ratio of actual/potential evap -!!\param[out] etp real, potential evaporation (\f$W/m^2\f$) +!!\param[out] drip real, through-fall of precip and/or dew in excess of canopy water-holding capacity (\f$m\f$) +!!\param[out] dew real, dewfall (or frostfall for t<273.15) (\f$m\f$) +!!\param[out] beta real, ratio of actual/potential evap +!!\param[out] etp real, potential evaporation (\f$W/m^2\f$) !!\param[out] ssoil real, upward soil heat flux (\f$W/m^2\f$) -!!\param[out] flx1 real, precip-snow sfc flux (\f$W/m^2\f$) -!!\param[out] flx2 real, freezing rain latent heat flux (\f$W/m^2\f$) -!!\param[out] flx3 real, phase-change heat flux from snowmelt (\f$W/m^2\f$) -!!\param[out] runoff1 real, surface runoff (\f$ms^{-1}\f$) not infiltrating sfc -!!\param[out] runoff2 real, sub sfc runoff (\f$ms^{-1}\f$) (baseflow) -!!\param[out] runoff3 real, excess of porosity for a given soil layer +!!\param[out] flx1 real, precip-snow sfc flux (\f$W/m^2\f$) +!!\param[out] flx2 real, freezing rain latent heat flux (\f$W/m^2\f$) +!!\param[out] flx3 real, phase-change heat flux from snowmelt (\f$W/m^2\f$) +!!\param[out] runoff1 real, surface runoff (\f$ms^{-1}\f$) not infiltrating sfc +!!\param[out] runoff2 real, sub sfc runoff (\f$ms^{-1}\f$) (baseflow) +!!\param[out] runoff3 real, excess of porosity for a given soil layer !!\param[out] snomlt real, snow melt (\f$m\f$) (water equivalent) !!\param[out] sncovr real, fractional snow cover -!!\param[out] rc real, canopy resistance (s/m) -!!\param[out] pc real, plant coeff (fraction) where pc*etp=transpi -!!\param[out] rsmin real, minimum canopy resistance (s/m) -!!\param[out] xlai real, leaf area index (dimensionless) -!!\param[out] rcs real, incoming solar rc factor (dimensionless) -!!\param[out] rct real, air temperature rc factor (dimensionless) -!!\param[out] rcq real, atoms vapor press deficit rc factor -!!\param[out] rcsoil real, soil moisture rc factor (dimensionless) -!!\param[out] soilw real, available soil moisture in root zone +!!\param[out] rc real, canopy resistance (s/m) +!!\param[out] pc real, plant coeff (fraction) where pc*etp=transpi +!!\param[out] rsmin real, minimum canopy resistance (s/m) +!!\param[out] xlai real, leaf area index (dimensionless) +!!\param[out] rcs real, incoming solar rc factor (dimensionless) +!!\param[out] rct real, air temperature rc factor (dimensionless) +!!\param[out] rcq real, atoms vapor press deficit rc factor +!!\param[out] rcsoil real, soil moisture rc factor (dimensionless) +!!\param[out] soilw real, available soil moisture in root zone !!\param[out] soilm real, total soil column moisture (frozen+unfrozen) (\f$m\f$) -!!\param[out] smcwlt real, wilting point (volumetric) -!!\param[out] smcdry real, dry soil moisture threshold (volumetric) -!!\param[out] smcref real, soil moisture threshold (volumetric) -!!\param[out] smcmax real, porosity (sat val of soil mois) +!!\param[out] smcwlt real, wilting point (volumetric) +!!\param[out] smcdry real, dry soil moisture threshold (volumetric) +!!\param[out] smcref real, soil moisture threshold (volumetric) +!!\param[out] smcmax real, porosity (sat val of soil mois) !!\section general_sflx GFS Noah LSM General Algorithm !! @{ ! subroutine sflx & ! --- inputs: @@ -272,7 +272,7 @@ subroutine gfssflx & ! --- inpu ! --- constant parameters: ! *** note: some of the constants are different in subprograms and need to ! be consolidated with the standard def in module physcons at sometime -! at the present time, those diverse values are kept temperately to +! at the present time, those diverse values are kept temperately to ! provide the same result as the original codes. -- y.t.h. may09 integer, parameter :: nsold = 4 ! max soil layers @@ -305,7 +305,7 @@ subroutine gfssflx & ! --- inpu real (kind=kind_phys), intent(in) :: ffrozp, dt, zlvl, lwdn, & & sldpth(nsoil), swdn, swnet, sfcems, sfcprs, sfctmp, & & sfcspd, prcp, q2, q2sat, dqsdt2, th2, shdmin, alb, snoalb, & - & bexpp, xlaip + & bexpp, xlaip ! --- input/outputs: real (kind=kind_phys), intent(inout) :: tbot, cmc, t1, sneqv, & @@ -354,7 +354,7 @@ subroutine gfssflx & ! --- inpu ! note - for open-sea, sflx should *not* have been called. set green ! vegetation fraction (shdfac) = 0. !> - For open-sea, sea-ice and glacial-ice cases, sflx() should not have -!! been called (set green vegetation fraction (shdfac) =0.) +!! been called (set green vegetation fraction (shdfac) =0.) ice = icein if(ivegsrc == 2) then @@ -385,7 +385,7 @@ subroutine gfssflx & ! --- inpu else -! --- ... calculate depth (negative) below ground from top skin sfc to +! --- ... calculate depth (negative) below ground from top skin sfc to ! bottom of each soil layer. ! note - sign of zsoil is negative (denoting below ground) @@ -395,8 +395,8 @@ subroutine gfssflx & ! --- inpu end do endif ! end if_ice_block - -! --- ... next is crucial call to set the land-surface parameters, + +! --- ... next is crucial call to set the land-surface parameters, ! including soil-type and veg-type dependent parameters. ! set shdfac=0.0 for bare soil surfaces @@ -478,8 +478,8 @@ subroutine gfssflx & ! --- inpu enddo endif -!> - If input snowpack (\a sneqv) is nonzero, then call csnow() to compute -!! snow density (\a sndens) and snow thermal conductivity (\a sncond). +!> - If input snowpack (\a sneqv) is nonzero, then call csnow() to compute +!! snow density (\a sndens) and snow thermal conductivity (\a sncond). ! (note that csnow is a function subroutine) if (sneqv .eq. 0.0) then @@ -499,9 +499,9 @@ subroutine gfssflx & ! --- inpu endif !> - Determine if it's precipitating and what kind of precipitation it is. -!! if it's precipitating and the air temperature is colder than \f$0^oC\f$, -!! it's snowing! if it's precipitating and the air temperature is warmer than -!! \f$0^oC\f$, but the ground temperature is colder than \f$0^oC\f$, freezing +!! if it's precipitating and the air temperature is colder than \f$0^oC\f$, +!! it's snowing! if it's precipitating and the air temperature is warmer than +!! \f$0^oC\f$, but the ground temperature is colder than \f$0^oC\f$, freezing !! rain is presumed to be falling. if (prcp > 0.0) then @@ -513,7 +513,7 @@ subroutine gfssflx & ! --- inpu endif !> - If either precipitation flag (\a snowng, \a frzgra) is set as true: -! determine new snowfall (converting precipitation rate from +! determine new snowfall (converting precipitation rate from ! \f$kg m^{-2} s^{-1}\f$ to a liquid equiv snow depth in meters) ! and add it to the existing snowpack. !> - Since all precip is added to snowpack, no precip infiltrates @@ -525,8 +525,8 @@ subroutine gfssflx & ! --- inpu sneqv = sneqv + sn_new prcp1 = 0.0 -!> - Call snow_new() to update snow density based on new snowfall, -!! using old and new snow. +!> - Call snow_new() to update snow density based on new snowfall, +!! using old and new snow. call snow_new ! --- inputs: ! ! ( sfctmp, sn_new, ! @@ -550,7 +550,7 @@ subroutine gfssflx & ! --- inpu endif ! end if_snowng_block -!> - Determine snowcover fraction and albedo fraction over sea-ice, +!> - Determine snowcover fraction and albedo fraction over sea-ice, !! glacial-ice, and land. For nonzero snow depth over land case: if (ice /= 0) then @@ -607,18 +607,18 @@ subroutine gfssflx & ! --- inpu ! --- ... next calculate the subsurface heat flux, which first requires ! calculation of the thermal diffusivity. treatment of the -! latter follows that on pages 148-149 from "heat transfer in -! cold climates", by v. j. lunardini (published in 1981 -! by van nostrand reinhold co.) i.e. treatment of two contiguous -! "plane parallel" mediums (namely here the first soil layer -! and the snowpack layer, if any). this diffusivity treatment -! behaves well for both zero and nonzero snowpack, including the +! latter follows that on pages 148-149 from "heat transfer in +! cold climates", by v. j. lunardini (published in 1981 +! by van nostrand reinhold co.) i.e. treatment of two contiguous +! "plane parallel" mediums (namely here the first soil layer +! and the snowpack layer, if any). this diffusivity treatment +! behaves well for both zero and nonzero snowpack, including the ! limit of very thin snowpack. this treatment also eliminates -! the need to impose an arbitrary upper bound on subsurface +! the need to impose an arbitrary upper bound on subsurface ! heat flux when the snowpack becomes extremely thin. ! --- ... first calculate thermal diffusivity of top soil layer, using -! both the frozen and liquid soil moisture, following the +! both the frozen and liquid soil moisture, following the ! soil thermal diffusivity function of peters-lidard et al. ! (1998,jas, vol 55, 1209-1224), which requires the specifying ! the quartz content of the given soil class (see routine redprm) @@ -636,15 +636,15 @@ subroutine gfssflx & ! --- inpu if ( vegtyp == 13 ) df1=3.24 endif -!> - Add subsurface heat flux reduction effect from the -!! overlying green canopy, adapted from section 2.1.2 of +!> - Add subsurface heat flux reduction effect from the +!! overlying green canopy, adapted from section 2.1.2 of !! \cite peters-lidard_et_al_1997. df1 = df1 * exp( sbeta*shdfac ) endif ! end if_ice_block -! --- ... finally "plane parallel" snowpack effect following +! --- ... finally "plane parallel" snowpack effect following ! v.j. linardini reference cited above. note that dtot is ! combined depth of snowdepth and thickness of first soil layer @@ -685,7 +685,7 @@ subroutine gfssflx & ! --- inpu endif ! end if_sneqv_block -!> - For uncoupled mode, call snowz0() to calculate surface roughness +!> - For uncoupled mode, call snowz0() to calculate surface roughness !! (\a z0) over snowpack using snow condition from the previous timestep. ! if (couple == 0) then ! uncoupled mode @@ -707,17 +707,17 @@ subroutine gfssflx & ! --- inpu ! --- ... next call routine sfcdif to calculate the sfc exchange coef (ch) ! for heat and moisture. -! note - comment out call sfcdif, if sfcdif already called in calling +! note - comment out call sfcdif, if sfcdif already called in calling ! program (such as in coupled atmospheric model). ! - do not call sfcdif until after above call to redprm, in case ! alternative values of roughness length (z0) and zilintinkevich ! coef (czil) are set there via namelist i/o. ! - routine sfcdif returns a ch that represents the wind spd times ! the "original" nondimensional "ch" typical in literature. hence -! the ch returned from sfcdif has units of m/s. the important +! the ch returned from sfcdif has units of m/s. the important ! companion coefficient of ch, carried here as "rch", is the ch ! from sfcdif times air density and parameter "cp". "rch" is -! computed in "call penman". rch rather than ch is the coeff +! computed in "call penman". rch rather than ch is the coeff ! usually invoked later in eqns. ! - sfcdif also returns the surface exchange coefficient for momentum, ! cm, also known as the surface drage coefficient, but cm is not @@ -783,7 +783,7 @@ subroutine gfssflx & ! --- inpu if (shdfac > 0.) then -! --- ... frozen ground extension: total soil water "smc" was replaced +! --- ... frozen ground extension: total soil water "smc" was replaced ! by unfrozen soil water "sh2o" in call to canres below call canres @@ -804,7 +804,7 @@ subroutine gfssflx & ! --- inpu if (sneqv .eq. 0.0) then !> - For no snowpack is present, call nopac() to calculate soil moisture !! and heat flux values and update soil moisture contant and soil heat -!! content values. +!! content values. call nopac ! --- inputs: ! ! ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, ! @@ -836,7 +836,7 @@ subroutine gfssflx & ! --- inpu ! ett, snomlt, drip, dew, flx1, flx3, esnow ) ! endif -!> - Noah LSM post-processing: +!> - Noah LSM post-processing: !> - Calculate sensible heat (h) for return to parent model. sheat = -(ch*cp1*sfcprs) / (rd1*t2v) * (th2 - t1) @@ -870,12 +870,12 @@ subroutine gfssflx & ! --- inpu !! - ssoil>0: warm the surface (night time) !! - ssoil<0: cool the surface (day time) - ssoil = -1.0 * ssoil + ssoil = -1.0 * ssoil if (ice == 0) then !> - For the case of land (but not glacial-ice): -!! convert runoff3 (internal layer runoff from supersat) from \f$m\f$ +!! convert runoff3 (internal layer runoff from supersat) from \f$m\f$ !! to \f$ms^-1\f$ and add to subsurface runoff/baseflow (runoff2). !! runoff2 is already a rate at this point. @@ -892,7 +892,7 @@ subroutine gfssflx & ! --- inpu runoff1 = snomlt / dt endif -!> - Calculate total column soil moisture in meters (soilm) and root-zone +!> - Calculate total column soil moisture in meters (soilm) and root-zone !! soil moisture availability (fraction) relative to porosity/saturation. soilm = -1.0 * smc(1) * zsoil(1) @@ -1228,7 +1228,7 @@ end subroutine csnow !----------------------------------- !>\ingroup Noah_LSM -!> This subroutine calculates soil moisture and heat flux values and +!> This subroutine calculates soil moisture and heat flux values and !! update soil moisture content and soil heat content values for the !! case when no snow pack is present. subroutine nopac @@ -1605,7 +1605,7 @@ end subroutine penman !----------------------------------- !>\ingroup Noah_LSM -!> This subroutine internally sets default values or optionally read-in +!> This subroutine internally sets default values or optionally read-in !! via namelist i/o, all soil and vegetation parateters requied for the execusion !! of the Noah LSM. subroutine redprm @@ -2257,7 +2257,7 @@ end subroutine snfrac !----------------------------------- !>\ingroup Noah_LSM -!> This subroutine calculates soil moisture and heat flux values and +!> This subroutine calculates soil moisture and heat flux values and !! update soil moisture content and soil heat content values for the !! case when a snow pack is present. subroutine snopac @@ -2566,7 +2566,9 @@ subroutine snopac ! so with snoexp = 2.0 (>1), surface skin temperature is higher than ! for the linear case (snoexp = 1). - t1 = tfreez * sncovr**snoexp + t12 * (1.0 - sncovr**snoexp) +! t1 = tfreez * sncovr**snoexp + t12 * (1.0 - sncovr**snoexp) + t1 = tfreez * max(0.01,sncovr**snoexp) + & + & t12 * (1.0 - max(0.01,sncovr**snoexp)) beta = 1.0 ssoil = df1 * (t1 - stc(1)) / dtot @@ -2620,7 +2622,7 @@ subroutine snopac ex = sneqv / dt flx3 = ex * 1000.0 * lsubf - snomlt = sneqv + snomlt = sneqv sneqv = 0.0 endif ! end if_sneqv-snomlt_block @@ -3577,8 +3579,8 @@ end subroutine smflx !----------------------------------- !>\ingroup Noah_LSM -!> This subroutine calculates compaction of a snowpack under conditions of -!! increasing snow density, as obtained from an approximate solution of +!> This subroutine calculates compaction of a snowpack under conditions of +!! increasing snow density, as obtained from an approximate solution of !! E. Anderson's differential equation (3.29),NOAA technical report NWS 19, !! by Victor Koren, 03/25/95. subroutine will return new values of \a snowh !! and \a sndens . @@ -3824,7 +3826,7 @@ end subroutine devap !! New version (June 2001): much faster and more accurate Newton iteration !! achieved by first taking log of eqn cited above -- less than 4 (typically !! 1 or 2) iterations achieves convergence. Also, explicit 1-step solution -!! option for special case of paramter ck=0, which reduces the orginal +!! option for special case of paramter ck=0, which reduces the orginal !! implicit equation to a simpler explicit form, known as the "flerchinger eqn". !! Improved handling of solution in the limit of freezing point temperature t0. subroutine frh2o & @@ -3979,7 +3981,7 @@ end subroutine frh2o !----------------------------------- !>\ingroup Noah_LSM -!> This subroutine calculates the right hand side of the time tendency +!> This subroutine calculates the right hand side of the time tendency !! term of the soil thermal diffusion equation. Also to compute (prepare) !! the matrix coefficients for the tri-diagonal matrix of the implicit time !! scheme. @@ -4854,7 +4856,7 @@ end subroutine snksrc !>\ingroup Noah_LSM !> This subroutine calculates the right hand side of the time tendency !! term of the soil water diffusion equation. Also to compute -!! (prepare) the matrix coefficients for the tri-diagonal matrix of +!! (prepare) the matrix coefficients for the tri-diagonal matrix of !! the implicit time scheme. subroutine srt & !................................... @@ -5316,7 +5318,7 @@ end subroutine sstep !----------------------------------- !>\ingroup Noah_LSM -!> This subroutine calculates temperature on the boundary of the +!> This subroutine calculates temperature on the boundary of the !! layer by interpolation of the middle layer temperatures. subroutine tbnd & !...................................