[Dart-dev] DART/branches Revision: 11462

[dart at ucar.edu]

thoar at ucar.edu
2017-04-12 18:46:54 -0600 (Wed, 12 Apr 2017)
393
this was called scripts_backup/model_mod_140418.f90  (2014 04 18)
original note:
    Add coefg_2003_10D into the original model_mod;
    Replace Ally's ss model with a radiative transfer model using the scheme of Njoku and Chan., 2006. 
    Add subroutine to read soil texture data from restart file;
    Constrain updated first layer soil moisture within 0-porosity (set to 0.5 if missing);




Modified: DART/branches/Tb/models/clm/model_mod.f90
===================================================================
--- DART/branches/Tb/models/clm/model_mod.f90	2017-04-13 00:43:51 UTC (rev 11461)
+++ DART/branches/Tb/models/clm/model_mod.f90	2017-04-13 00:46:54 UTC (rev 11462)
@@ -116,7 +116,7 @@
 character(len=32 ), parameter :: revision = "$Revision$"
 character(len=128), parameter :: revdate  = "$Date$"
 
-character(len=256) :: string1, string2, string3
+character(len=256) :: string1, string2, string3, string4
 logical, save :: module_initialized = .false.
 
 ! Storage for a random sequence for perturbing a single initial state
@@ -158,11 +158,13 @@
 character(len=256) :: clm_restart_filename = 'clm_restart.nc'
 character(len=256) :: clm_history_filename = 'clm_history.nc'
 character(len=256) :: casename = 'clm_dart'
+character(len=256) :: coefg_nc = 'coefg_amsre2003_10D.nc'
 
 character(len=obstypelength) :: clm_state_variables(max_state_variables*num_state_table_columns) = ' '
 
 namelist /model_nml/            &
-   casename,                    &
+   casename,                    & 
+   coefg_nc,                    &
    clm_restart_filename,        &
    clm_history_filename,        &
    output_state_vector,         &
@@ -203,6 +205,20 @@
 ! Properties required for a snow column
 !----------------------------------------------------------------------
 
+
+type soilprops
+   private
+   integer  :: nlayers  ! snow layers
+   real(r4) :: ssm      ! aux_ins(1) surface soil moisture [m3m-3]
+   real(r4) :: stg      ! aux_ins(2) surface layer temperature [k]   
+   real(r4) :: sat      ! aux_ins(3) topsoil porosity
+   real(r4) :: ssand    ! aux_ins(4) soil sand content percentage (0~100)
+   real(r4) :: sclay    ! aux_ins(5) soil clay content percentage (0~100)
+   real(r4) :: scoefg   ! aux_ins(6) background value for coefficient "g"
+end type soilprops
+
+type(soilprops) :: soilcolumn
+
 type snowprops
    private
    integer  :: nlayers         ! aux_ins(1)
@@ -2328,6 +2344,7 @@
 
 ! BOMBPROOFING - check for a vertical dimension for this variable
 if (progvar(ivar)%maxlevels > 1) then
+   write(*,*)'progvar(ivar)%maxlevels = ',progvar(ivar)%maxlevels  !====Long
    write(string1, *)'Variable '//trim(varstring)//' cannot use this routine.'
    write(string2, *)'use get_grid_vertval() instead.'
    call error_handler(E_ERR,'compute_gridcell_value', string1, &
@@ -2732,7 +2749,8 @@
 integer, optional,        intent(in)  :: ncid
 
 integer :: i,j,ii, VarID
-real(r8), allocatable, dimension(:,:) :: org_array
+real(r8), allocatable, dimension(:,:) :: org_array, org_porosity
+real(r8), allocatable, dimension(:)   :: org_watsat
 
 if ( .not. module_initialized ) call static_init_model
 
@@ -2761,6 +2779,8 @@
 if (present(ncid)) then
 
    allocate(org_array(size(data_2d_array,1),size(data_2d_array,2)))
+   allocate(org_porosity(15,size(data_2d_array,2)))      !============Long
+   allocate(org_watsat(size(data_2d_array,2)))           !============Long
 
    call nc_check(nf90_inq_varid(ncid, progvar(ivar)%varname, VarID), &
             'vector_to_2d_prog_var', 'inq_varid '//trim(progvar(ivar)%varname))
@@ -2768,6 +2788,16 @@
    call nc_check(nf90_get_var(ncid, VarID, org_array), &
             'vector_to_2d_prog_var', 'get_var '//trim(progvar(ivar)%varname))
 
+   call nc_check(nf90_inq_varid(ncid, 'WATSAT', VarID), &
+            'vector_to_2d_prog_var', 'inq_varid WATSAT')
+
+   call nc_check(nf90_get_var(ncid, VarID, org_porosity), &
+            'vector_to_2d_prog_var', 'get_var WATSAT')
+
+   ! to convert top layer soil porosity with unit: m3/m3 to kg/m2
+   where((org_porosity>1)) org_porosity=0.5_r8
+   org_watsat=org_porosity(1,:)*(LEVGRND(6)+LEVGRND(7))/2*1000
+
    ! restoring the indeterminate original values
 
    where(data_2d_array == MISSING_R8 ) data_2d_array = org_array
@@ -2781,6 +2811,20 @@
    elseif (trim(progvar(ivar)%varname) == 'ZISNO') then
       where((data_2d_array < 0.0_r8)) data_2d_array = org_array
    elseif (trim(progvar(ivar)%varname) == 'H2OSOI_LIQ') then
+      !===========================================================Long
+      ! Currently, only update the first layer of soil, 
+      ! i.e., the 6th layer in levtot.
+      ! So replace other layers' value with original value.