! DART software - Copyright 2004 - 2013 UCAR. This open source software is ! provided by UCAR, "as is", without charge, subject to all terms of use at ! http://www.image.ucar.edu/DAReS/DART/DART_download ! ! $Id: model_mod.f90 7202 2014-10-08 15:29:19Z thoar $ module model_mod ! This is the interface between the Community Land Model (CLM) and DART. ! The following is the hierarchy as I see it: ! top ... a gridcell has one or more more landunits ! ... ... landunits have one or more columns ! ... ... ... columns may have one or more layers (snow, for instance) ! ... ... ... columns have one or more PFTs ! ... ... ... ... some PFTs have layers (radiance bands, for instance) ! ! landunit types ! 1 soil (natural vegation/bare ground) ! 2 glacier ! 3 lake ! 4 not used (shallow lake) ! 5 wetland ! 6 urban ! 7 ice (new glacier model) ! 8 crop (if using crop model) ! Modules that are absolutely required for use are listed use types_mod, only : r4, r8, SECPERDAY, MISSING_R8, & MISSING_I, MISSING_R4, rad2deg, deg2rad, PI, & obstypelength use time_manager_mod, only : time_type, set_time, get_time, set_date, get_date,& print_time, print_date, set_calendar_type, & operator(*), operator(+), operator(-), & operator(>), operator(<), operator(/), & operator(/=), operator(<=), operator(==) use location_mod, only : location_type, get_dist, query_location, & get_close_maxdist_init, get_close_type, & set_location, get_location, horiz_dist_only, & vert_is_undef, VERTISUNDEF, & vert_is_surface, VERTISSURFACE, & vert_is_level, VERTISLEVEL, & vert_is_pressure, VERTISPRESSURE, & vert_is_height, VERTISHEIGHT, & get_close_obs_init, get_close_obs, LocationDims use utilities_mod, only : register_module, error_handler, & E_ERR, E_WARN, E_MSG, logfileunit, get_unit, & nc_check, do_output, to_upper, & find_namelist_in_file, check_namelist_read, & file_exist, find_textfile_dims, file_to_text, & open_file, close_file use obs_kind_mod, only : KIND_SOIL_TEMPERATURE, & KIND_SOIL_MOISTURE, & KIND_LIQUID_WATER, & KIND_ICE, & KIND_SNOWCOVER_FRAC, & KIND_SNOW_THICKNESS, & KIND_LEAF_CARBON, & KIND_WATER_TABLE_DEPTH, & KIND_GEOPOTENTIAL_HEIGHT,& KIND_BRIGHTNESS_TEMPERATURE, & paramname_length, & get_raw_obs_kind_index use mpi_utilities_mod, only: my_task_id use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian use typesizes use netcdf implicit none private ! these routines must be public and you cannot change ! the arguments - they will be called *from* the DART code. public :: get_model_size, & adv_1step, & get_state_meta_data, & model_interpolate, & get_model_time_step, & static_init_model, & end_model, & init_time, & init_conditions, & nc_write_model_atts, & nc_write_model_vars, & pert_model_state, & get_close_maxdist_init, & get_close_obs_init, & get_close_obs, & ens_mean_for_model ! generally useful routines for various support purposes. ! the interfaces here can be changed as appropriate. public :: clm_to_dart_state_vector, & sv_to_restart_file, & get_clm_restart_filename, & get_state_time, & get_grid_vertval, & compute_gridcell_value, & gridcell_components, & DART_get_var, & get_model_time, & get_ncols_in_gridcell, & get_colids_in_gridcell ! version controlled file description for error handling, do not edit character(len=256), parameter :: source = & "$URL: https://proxy.subversion.ucar.edu/DAReS/DART/branches/Tb/models/clm/model_mod.f90 $" character(len=32 ), parameter :: revision = "$Revision: 7202 $" character(len=128), parameter :: revdate = "$Date: 2014-10-08 10:29:19 -0500 (Wed, 08 Oct 2014) $" character(len=256) :: string1, string2, string3 logical, save :: module_initialized = .false. ! Storage for a random sequence for perturbing a single initial state type(random_seq_type) :: random_seq !------------------------------------------------------------------ ! ! The DART state vector may consist of things like: ! ! SNOWDP aka "snow depth" ! frac_sno aka "snow cover fraction" ! leafc aka "leaf carbon" ! T_SOISNO aka "temperature of soil & snow" ! H2OSOI_LIQ aka "liquid water in soil & snow" ! H2OSOI_ICE aka "water equivalent of ice in soil & snow" ! DZSNO aka "snow layer thickness" ! ! The variables in the clm restart file that are used to create the ! DART state vector are specified in the input.nml:model_nml namelist. ! !------------------------------------------------------------------ integer, parameter :: LAKE = 3 ! Codes for restricting the range of a variable integer, parameter :: BOUNDED_NONE = 0 ! ... unlimited range integer, parameter :: BOUNDED_BELOW = 1 ! ... minimum, but no maximum integer, parameter :: BOUNDED_ABOVE = 2 ! ... maximum, but no minimum integer, parameter :: BOUNDED_BOTH = 3 ! ... minimum and maximum integer :: nfields integer, parameter :: max_state_variables = 40 integer, parameter :: num_state_table_columns = 6 character(len=obstypelength) :: variable_table(max_state_variables, num_state_table_columns) ! Codes for interpreting the columns of the variable_table integer, parameter :: VT_VARNAMEINDX = 1 ! ... variable name integer, parameter :: VT_KINDINDX = 2 ! ... DART kind integer, parameter :: VT_MINVALINDX = 3 ! ... minimum value if any integer, parameter :: VT_MAXVALINDX = 4 ! ... maximum value if any integer, parameter :: VT_ORIGININDX = 5 ! ... file of origin integer, parameter :: VT_STATEINDX = 6 ! ... update (state) or not ! things which can/should be in the model_nml integer :: assimilation_period_days = 0 integer :: assimilation_period_seconds = 60 real(r8) :: model_perturbation_amplitude = 0.2 logical :: output_state_vector = .true. integer :: debug = 0 ! turn up for more and more debug messages character(len=32) :: calendar = 'Gregorian' character(len=256) :: clm_restart_filename = 'clm_restart.nc' character(len=256) :: clm_history_filename = 'clm_history.nc' character(len=256) :: clm_vector_history_filename = 'clm_vector_history.nc' character(len=256) :: casename = 'clm_dart' character(len=256) :: lai_nc = '../rtm_lai_modis.nc' !kyh01192016 character(len=obstypelength) :: clm_variables(max_state_variables*num_state_table_columns) = ' ' real(r8),dimension(288,192,40) :: tran_can_v_grd !kyh03092016 (lon=288, lat=192, ensemble_size=40) real(r8),dimension(288,192,40) :: tran_can_h_grd !kyh03092016 (lon=288, lat=192, ensemble_size=40) namelist /model_nml/ & casename, & clm_restart_filename, & clm_history_filename, & clm_vector_history_filename, & output_state_vector, & assimilation_period_days, & ! for now, this is the timestep assimilation_period_seconds, & model_perturbation_amplitude,& calendar, & debug, & clm_variables, & lai_nc !kyh01192016 ! Everything needed to describe a variable type progvartype private character(len=NF90_MAX_NAME) :: varname character(len=NF90_MAX_NAME) :: long_name character(len=NF90_MAX_NAME) :: units character(len=obstypelength), dimension(NF90_MAX_VAR_DIMS) :: dimnames integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens integer :: numdims integer :: maxlevels integer :: xtype integer :: varsize ! prod(dimlens(1:numdims)) integer :: index1 ! location in dart state vector of first occurrence integer :: indexN ! location in dart state vector of last occurrence integer :: dart_kind integer :: rangeRestricted real(r8) :: minvalue real(r8) :: maxvalue integer :: spvalINT, missingINT real(r4) :: spvalR4, missingR4 real(r8) :: spvalR8, missingR8 logical :: has_fill_value ! intended for future use logical :: has_missing_value ! intended for future use character(len=paramname_length) :: kind_string character(len=512) :: origin ! the file it came from logical :: update end type progvartype type(progvartype), dimension(max_state_variables) :: progvar !---------------------------------------------------------------------- ! Properties required for a snow column !---------------------------------------------------------------------- type snowprops private integer :: nlayers ! aux_ins(1) real(r4) :: t_grnd ! aux_ins(2) ground temperature [K] real(r4) :: soilsat ! aux_ins(3) soil saturation [fraction] real(r4) :: soilporos ! aux_ins(4) soil porosity [fraction] real(r4) :: propconst ! aux_ins(5) proportionality between grain size & correlation length. !----------------------------------------------------------------------------------------------kyh04032014 real(r4) :: soil_liq ! aux_ins(6) soil liquid water content (fraction) real(r4) :: soil_ice ! aux_ins(7) soil ice contet (fraction) real(r4) :: sandf ! aux_ins(8) sand fraction (fraction) real(r4) :: clayf ! aux_ins(9) clay fraction (fraction) real(r4) :: Tcanopy ! aux_ins(10) vegetation physical temperature (K) real(r4) :: lai_vege ! aux_ins(11) leaf area index real(r4) :: forc_pbot ! aux_ins(12) ground level pressure (Pa) real(r4) :: forc_t ! aux_ins(13) ground level air temperature (K) real(r4) :: forc_rh ! aux_ins(14) relative humidity real(r4) :: vf ! aux_ins(15) vegetated area fraction !----------------------------------------------------------------------------------------------kyh04032014 !----------------------------------------------------------------------------------------------RTM parameters !kyh11202014 real(r4) :: stickiness !snowpack stickiness (-) real(r4) :: b_prime !vegetation RTM coefficient (-) real(r4) :: x_lambda !vegetation RTM coefficient (-) !----------------------------------------------------------------------------------------------RTM parameters !kyh11202014 real(r4) :: scatt_albedo !vegetation single scattering albedo (-) !kyh03052016 integer :: nprops ! [thickness, density, diameter, liqwater, temperature] real(r4), pointer, dimension(:) :: thickness ! LAYER THICKNESS [M] real(r4), pointer, dimension(:) :: density ! LAYER DENSITY [KG/M3] real(r4), pointer, dimension(:) :: grain_radius ! LAYER GRAIN RADIUS [M] real(r4), pointer, dimension(:) :: liquid_water ! LAYER LIQUID WATER CONTENT [FRAC] real(r4), pointer, dimension(:) :: temperature ! LAYER TEMPERATURE [K] end type snowprops type(snowprops) :: snowcolumn !---------------------------------------------------------------------- ! how many and which columns are in each gridcell !---------------------------------------------------------------------- type gridcellcolumns ! given a gridcell, which columns contribute private integer :: ncols integer, pointer, dimension(:) :: columnids integer :: npfts integer, pointer, dimension(:) :: pftids end type gridcellcolumns type(gridcellcolumns), allocatable, dimension(:,:), target :: gridCellInfo !------------------------------------------------------------------------------ ! Things that come from the CLM history file. ! ! The LON,LAT arrays store the longitude and latitude of grid cell centers. ! For the FV cores, there actually _are_ cells CENTERED at the poles. integer :: nlon = -1 integer :: nlonatm = -1 integer :: nlonrof = -1 integer :: nlat = -1 integer :: nlatatm = -1 integer :: nlatrof = -1 integer :: nlevgrnd = -1 ! Number of soil levels real(r8), allocatable :: LEVGRND(:) real(r8), allocatable :: LON(:) real(r8), allocatable :: LAT(:) real(r8), allocatable :: AREA1D(:), LANDFRAC1D(:) ! unstructured grid real(r8), allocatable :: AREA2D(:,:), LANDFRAC2D(:,:) ! 2D grid logical :: unstructured = .false. !------------------------------------------------------------------------------ ! Things that come from the CLM restart file. ! ! These are the 'sparse' type arrays pertaining to the land gridcells. ! ! ZSNO contains the height of the middle of each snow layer; ! ZISNO contains the height of the top of each snow layer; ! DZSNO tells the snow thickness of each layer. ! snow heights are stored as negative values ! Unlike the soil layers, snow layer thickness, as well as snow layer depth, ! may change as time goes on (due to snow metamorphism, overburden and the like). ! So there's no uniform levsno as levgrnd coordinate variable. ! The packing order is: Z-LON-LAT, Z moving fastest. ! all levels at a location, then ! scan along longitudes, then ! move to next latitude. integer :: Ngridcell = -1 ! Number of gridcells containing land integer :: Nlandunit = -1 ! Number of land units integer :: Ncolumn = -1 ! Number of columns integer :: Npft = -1 ! Number of plant functional types integer :: Nlevlak = -1 ! Number of integer :: Nlevsno = -1 ! Number of snow levels integer :: Nlevsno1 = -1 ! Number of snow level ... interfaces? integer :: Nlevtot = -1 ! Number of total levels integer :: Nnumrad = -1 ! Number of integer :: Nrtmlon = -1 ! Number of river transport model longitudes integer :: Nrtmlat = -1 ! Number of river transport model latitudes integer, allocatable, dimension(:) :: grid1d_ixy, grid1d_jxy ! 2D lon/lat index of corresponding gridcell integer, allocatable, dimension(:) :: land1d_ixy, land1d_jxy ! 2D lon/lat index of corresponding gridcell integer, allocatable, dimension(:) :: cols1d_ixy, cols1d_jxy ! 2D lon/lat index of corresponding gridcell integer, allocatable, dimension(:) :: pfts1d_ixy, pfts1d_jxy ! 2D lon/lat index of corresponding gridcell real(r8), allocatable, dimension(:) :: land1d_wtxy ! landunit weight relative to corresponding gridcell real(r8), allocatable, dimension(:) :: cols1d_wtxy ! column weight relative to corresponding gridcell real(r8), allocatable, dimension(:) :: pfts1d_wtxy ! pft weight relative to corresponding gridcell real(r8), allocatable, dimension(:) :: levtot real(r8), allocatable, dimension(:,:):: zsno ! (column,levsno) ... snow layer midpoint integer, allocatable, dimension(:) :: cols1d_ityplun ! columntype ... lake, forest, city ... !------------------------------------------------------------------------------ ! These are the metadata arrays that are the same size as the state vector. real(r8), allocatable, dimension(:) :: ens_mean ! may be needed for forward ops integer, allocatable, dimension(:) :: lonixy ! longitude index of parent gridcell integer, allocatable, dimension(:) :: latjxy ! latitude index of parent gridcell real(r8), allocatable, dimension(:) :: levels ! depth real(r8), allocatable, dimension(:) :: landarea ! land area ... 'support' ... 'weight' !------------------------------------------------------------------------------ ! set this to true if you want to print out the current time ! after each N observations are processed, for benchmarking. logical :: print_timestamps = .false. integer :: print_every_Nth = 10000 !------------------------------------------------------------------ ! module storage !------------------------------------------------------------------ integer :: model_size ! the state vector length type(time_type) :: model_time ! valid time of the model state type(time_type) :: model_timestep ! smallest time to adv model INTERFACE vector_to_prog_var MODULE PROCEDURE vector_to_1d_prog_var MODULE PROCEDURE vector_to_2d_prog_var END INTERFACE INTERFACE DART_get_var MODULE PROCEDURE get_var_1d MODULE PROCEDURE get_var_2d MODULE PROCEDURE get_var_3d END INTERFACE INTERFACE get_state_time MODULE PROCEDURE get_state_time_ncid MODULE PROCEDURE get_state_time_fname END INTERFACE contains !================================================================== ! All the REQUIRED interfaces come first - just by convention. !================================================================== function get_model_size() !------------------------------------------------------------------ ! Returns the size of the model as an integer. ! Required for all applications. integer :: get_model_size if ( .not. module_initialized ) call static_init_model get_model_size = model_size end function get_model_size subroutine adv_1step(x, time) !------------------------------------------------------------------ ! CLM is never advanced by DART. ! If we get here, something is wrong and we should stop right away. real(r8), intent(inout) :: x(:) type(time_type), intent(in) :: time if ( .not. module_initialized ) call static_init_model if (do_output()) then call print_time(time,'NULL interface adv_1step (no advance) DART time is') call print_time(time,'NULL interface adv_1step (no advance) DART time is',logfileunit) endif write(string1,*)'DART should not be trying to advance CLM' call error_handler(E_ERR,'adv_1step',string1,source,revision,revdate) ! just so suppress compiler warnings. code unreachable x(:) = MISSING_R8 end subroutine adv_1step subroutine get_state_meta_data(indx, location, var_type) !------------------------------------------------------------------ ! Given an integer index into the state vector structure, returns the ! associated array indices for lat, lon, and height, as well as the type. ! Passed variables integer, intent(in) :: indx type(location_type) :: location integer, OPTIONAL, intent(out) :: var_type ! Local variables integer :: n ! Module variables ! LON ! LAT ! lonixy ! latjxy ! levels ! progvar if ( .not. module_initialized ) call static_init_model location = set_location( LON(lonixy(indx)), LAT(latjxy(indx)), levels(indx), VERTISHEIGHT) if (present(var_type)) then var_type = MISSING_I FINDTYPE : do n = 1,nfields if((indx >= progvar(n)%index1) .and. & (indx <= progvar(n)%indexN) ) then var_type = progvar(n)%dart_kind exit FINDTYPE endif enddo FINDTYPE if( var_type == MISSING_I ) then write(string1,*) 'Problem, cannot find base_offset, indx is: ', indx call error_handler(E_ERR,'get_state_meta_data',string1,source,revision,revdate) endif endif return end subroutine get_state_meta_data function get_model_time_step() !------------------------------------------------------------------ ! ! Returns the the time step of the model; the smallest increment ! in time that the model is capable of advancing the state in a given ! implementation. This interface is required for all applications. type(time_type) :: get_model_time_step if ( .not. module_initialized ) call static_init_model get_model_time_step = model_timestep end function get_model_time_step subroutine static_init_model() !------------------------------------------------------------------ ! ! Called to do one time initialization of the model. ! ! All the grid information comes from the initialization of ! the dart_clm_mod module. ! Local variables - all the important ones have module scope integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs character(len=obstypelength) :: dimname integer :: ncid, TimeDimID, VarID, dimlen, varsize integer :: iunit, io, ivar integer :: i, j, xi, xj, index1, indexN, indx integer :: ss, dd integer :: spvalINT real(r4) :: spvalR4 real(r8) :: spvalR8 if ( module_initialized ) return ! only need to do this once. ! Print module information to log file and stdout. call register_module(source, revision, revdate) ! Since this routine calls other routines that could call this routine ! we'll say we've been initialized pretty dang early. module_initialized = .true. ! Read the DART namelist for this model call find_namelist_in_file('input.nml', 'model_nml', iunit) read(iunit, nml = model_nml, iostat = io) call check_namelist_read(iunit, io, 'model_nml') ! Record the namelist values used for the run if (do_output()) call error_handler(E_MSG,'static_init_model','model_nml values are') if (do_output()) write(logfileunit, nml=model_nml) if (do_output()) write( * , nml=model_nml) !--------------------------------------------------------------- ! Set the time step ... causes clm namelists to be read. ! Ensures model_timestep is multiple of 'dynamics_timestep' call set_calendar_type( calendar ) ! comes from model_mod_nml model_time = get_state_time(clm_restart_filename) model_timestep = set_model_time_step() call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400) write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds' call error_handler(E_MSG,'static_init_model',string1) !--------------------------------------------------------------- ! The CLM history file (h0?) has the 'superset' of information. ! The CLM restart files are intentionally lean and, in so doing, ! do not have the full lat/lon arrays nor any depth information. ! ncid = 0; ! signal that netcdf file is closed call get_history_dims(ncid, clm_history_filename, 'open', nlon, nlat, nlevgrnd ) if (unstructured) then allocate( AREA1D(nlon), LANDFRAC1D(nlon) ) else allocate( AREA2D(nlon,nlat), LANDFRAC2D(nlon,nlat) ) endif allocate(LON(nlon), LAT(nlat), LEVGRND(nlevgrnd)) call get_full_grid(ncid, clm_history_filename, 'close') ncid = 0; ! signal that netcdf file is closed !--------------------------------------------------------------- ! The CLM grid in a restart file is fundamentally a sparse matrix ! representation that lacks the native grid dimensions. ! The full lat/lon arrays do not exist in the restart files. ! only the grid cells that contain land are preserved. call get_sparse_dims(ncid, clm_restart_filename, 'open') allocate(grid1d_ixy(Ngridcell), grid1d_jxy(Ngridcell)) allocate(land1d_ixy(Nlandunit), land1d_jxy(Nlandunit), land1d_wtxy(Nlandunit)) allocate(cols1d_ixy(Ncolumn), cols1d_jxy(Ncolumn)) allocate(cols1d_wtxy(Ncolumn), cols1d_ityplun(Ncolumn)) allocate(pfts1d_ixy(Npft), pfts1d_jxy(Npft) , pfts1d_wtxy(Npft)) allocate(levtot(Nlevtot)) if (Nlevsno > 0) allocate(zsno(Nlevsno,Ncolumn)) call get_sparse_geog(ncid, clm_restart_filename, 'close') !--------------------------------------------------------------- ! Generate list of columns in each gridcell allocate(gridCellInfo(nlon,nlat)) call SetLocatorArrays() !--------------------------------------------------------------- ! Compile the list of clm variables to use in the creation ! of the DART state vector. call parse_variable_table( clm_variables, nfields, variable_table ) ! Compute the offsets into the state vector for the start of each ! variable type. Requires reading shapes from the clm restart file. ! Record the extent of the data type in the state vector. index1 = 1; indexN = 0; do ivar = 1, nfields ! convey the information in the variable_table to each progvar. call SetVariableAttributes(ivar) ! Open the file for each variable and get dimensions, etc. call nc_check(nf90_open(trim(progvar(ivar)%origin), NF90_NOWRITE, ncid), & 'static_init_model','open '//trim(progvar(ivar)%origin)) ! File is not required to have a time dimension io = nf90_inq_dimid(ncid, 'time', TimeDimID) if (io /= NF90_NOERR) TimeDimID = MISSING_I string2 = trim(progvar(ivar)%origin)//' '//trim(progvar(ivar)%varname) call nc_check(nf90_inq_varid(ncid, trim(progvar(ivar)%varname), VarID), & 'static_init_model', 'inq_varid '//trim(string2)) call nc_check(nf90_inquire_variable(ncid, VarID, dimids=dimIDs, & ndims=progvar(ivar)%numdims, xtype=progvar(ivar)%xtype), & 'static_init_model', 'inquire '//trim(string2)) ! If the long_name and/or units attributes are set, get them. ! They are not REQUIRED to exist but are nice to use if they are present. if( nf90_inquire_attribute( ncid, VarID, 'long_name') == NF90_NOERR ) then call nc_check( nf90_get_att(ncid, VarID, 'long_name' , progvar(ivar)%long_name), & 'static_init_model', 'get_att long_name '//trim(string2)) else progvar(ivar)%long_name = progvar(ivar)%varname endif if( nf90_inquire_attribute( ncid, VarID, 'units') == NF90_NOERR ) then call nc_check( nf90_get_att(ncid, VarID, 'units' , progvar(ivar)%units), & 'static_init_model', 'get_att units '//trim(string2)) else progvar(ivar)%units = 'unknown' endif ! Saving any FillValue, missing_value attributes so I can use it when I read and write ... ! CESM1_1_1 ... no attributes in the restart file for rank1 or greater ! variables. if (progvar(ivar)%xtype == NF90_INT) then if (nf90_get_att(ncid, VarID, '_FillValue' , spvalINT) == NF90_NOERR) then progvar(ivar)%spvalINT = spvalINT progvar(ivar)%has_fill_value = .true. endif if (nf90_get_att(ncid, VarID, 'missing_value' , spvalINT) == NF90_NOERR) then progvar(ivar)%missingINT = spvalINT progvar(ivar)%has_missing_value = .true. endif elseif (progvar(ivar)%xtype == NF90_FLOAT) then if (nf90_get_att(ncid, VarID, '_FillValue' , spvalR4) == NF90_NOERR) then progvar(ivar)%spvalR4 = spvalR4 progvar(ivar)%has_fill_value = .true. endif if (nf90_get_att(ncid, VarID, 'missing_value' , spvalR4) == NF90_NOERR) then progvar(ivar)%missingR4 = spvalR4 progvar(ivar)%has_missing_value = .true. endif elseif (progvar(ivar)%xtype == NF90_DOUBLE) then if (nf90_get_att(ncid, VarID, '_FillValue' , spvalR8) == NF90_NOERR) then progvar(ivar)%spvalR8 = spvalR8 progvar(ivar)%has_fill_value = .true. endif if (nf90_get_att(ncid, VarID, 'missing_value' , spvalR8) == NF90_NOERR) then progvar(ivar)%missingR8 = spvalR8 progvar(ivar)%has_missing_value = .true. endif endif varsize = 1 dimlen = 1 DimensionLoop : do i = 1,progvar(ivar)%numdims write(string1,'(''inquire dimension'',i2,A)') i,trim(string2) call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), name=dimname, len=dimlen), & 'static_init_model', string1) ! Only reserve space for a single time slice if (dimIDs(i) == TimeDimID) dimlen = 1 progvar(ivar)%dimlens( i) = dimlen progvar(ivar)%dimnames(i) = dimname varsize = varsize * dimlen enddo DimensionLoop progvar(ivar)%varsize = varsize progvar(ivar)%index1 = index1 progvar(ivar)%indexN = index1 + varsize - 1 index1 = index1 + varsize ! sets up for next variable if ((debug > 0) .and. do_output()) then write(logfileunit,*) write(logfileunit,*) trim(progvar(ivar)%varname),' variable number ',ivar write(logfileunit,*) ' filename ',trim(progvar(ivar)%origin) write(logfileunit,*) ' update ',progvar(ivar)%update write(logfileunit,*) ' long_name ',trim(progvar(ivar)%long_name) write(logfileunit,*) ' units ',trim(progvar(ivar)%units) write(logfileunit,*) ' xtype ',progvar(ivar)%xtype write(logfileunit,*) ' dimnames ',(/ (trim(progvar(ivar)%dimnames(i))//' ', i=1,progvar(ivar)%numdims ) /) write(logfileunit,*) ' dimlens ',progvar(ivar)%dimlens( 1:progvar(ivar)%numdims) write(logfileunit,*) ' numdims ',progvar(ivar)%numdims write(logfileunit,*) ' varsize ',progvar(ivar)%varsize write(logfileunit,*) ' index1 ',progvar(ivar)%index1 write(logfileunit,*) ' indexN ',progvar(ivar)%indexN write(logfileunit,*) ' dart_kind ',progvar(ivar)%dart_kind write(logfileunit,*) ' kind_string ',progvar(ivar)%kind_string write(logfileunit,*) ' spvalINT ',progvar(ivar)%spvalINT write(logfileunit,*) ' spvalR4 ',progvar(ivar)%spvalR4 write(logfileunit,*) ' spvalR8 ',progvar(ivar)%spvalR8 write(logfileunit,*) ' missingINT ',progvar(ivar)%missingINT write(logfileunit,*) ' missingR4 ',progvar(ivar)%missingR4 write(logfileunit,*) ' missingR8 ',progvar(ivar)%missingR8 write(logfileunit,*) ' has_fill_value ',progvar(ivar)%has_fill_value write(logfileunit,*) ' has_missing_value ',progvar(ivar)%has_missing_value write(logfileunit,*)' rangeRestricted ',progvar(ivar)%rangeRestricted write(logfileunit,*)' minvalue ',progvar(ivar)%minvalue write(logfileunit,*)' maxvalue ',progvar(ivar)%maxvalue write( * ,*) write( * ,*) trim(progvar(ivar)%varname),' variable number ',ivar write( * ,*) ' filename ',trim(progvar(ivar)%origin) write( * ,*) ' update ',progvar(ivar)%update write( * ,*) ' long_name ',trim(progvar(ivar)%long_name) write( * ,*) ' units ',trim(progvar(ivar)%units) write( * ,*) ' xtype ',progvar(ivar)%xtype write( * ,*) ' dimnames ',(/ (trim(progvar(ivar)%dimnames(i))//' ', i=1,progvar(ivar)%numdims ) /) write( * ,*) ' dimlens ',progvar(ivar)%dimlens( 1:progvar(ivar)%numdims) write( * ,*) ' numdims ',progvar(ivar)%numdims write( * ,*) ' varsize ',progvar(ivar)%varsize write( * ,*) ' index1 ',progvar(ivar)%index1 write( * ,*) ' indexN ',progvar(ivar)%indexN write( * ,*) ' dart_kind ',progvar(ivar)%dart_kind write( * ,*) ' kind_string ',progvar(ivar)%kind_string write( * ,*) ' spvalINT ',progvar(ivar)%spvalINT write( * ,*) ' spvalR4 ',progvar(ivar)%spvalR4 write( * ,*) ' spvalR8 ',progvar(ivar)%spvalR8 write( * ,*) ' missingINT ',progvar(ivar)%missingINT write( * ,*) ' missingR4 ',progvar(ivar)%missingR4 write( * ,*) ' missingR8 ',progvar(ivar)%missingR8 write( * ,*) ' has_fill_value ',progvar(ivar)%has_fill_value write( * ,*) ' has_missing_value ',progvar(ivar)%has_missing_value write( * ,*)' rangeRestricted ',progvar(ivar)%rangeRestricted write( * ,*)' minvalue ',progvar(ivar)%minvalue write( * ,*)' maxvalue ',progvar(ivar)%maxvalue endif call nc_check(nf90_close(ncid),'static_init_model','close '//trim(string2)) ncid = 0 enddo model_size = progvar(nfields)%indexN if ((debug > 0) .and. do_output()) then write(logfileunit, *) write(logfileunit,'("grid: nlon, nlat, nz =",2(1x,i6))') nlon, nlat write(logfileunit, *)'model_size = ', model_size write( * , *) write( * ,'("grid: nlon, nlat, nz =",2(1x,i6))') nlon, nlat write( * , *)'model_size = ', model_size endif allocate(ens_mean(model_size)) !--------------------------------------------------------------- ! Create the metadata arrays that are the same shape as the state vector. ! The metadata arrays will provide the ability to determine what grid cell is the parent ! of the state vector index in question ... as well as the actual surface area. ! This MUST stride through the state vector the same way the state vector is filled. allocate(lonixy(model_size), latjxy(model_size), levels(model_size), landarea(model_size)) ! Initialize all levels to surface. If there is a level, we will explicitly specify it. levels(:) = 0.0_r8 do ivar=1, nfields ! All variables are at the surface until proven otherwise. progvar(ivar)%maxlevels = 1 indx = progvar(ivar)%index1 ! 1D variables are usually from the restart file ! FIXME this same logic is used for 2D variables from the 'vector' file which ! has a singleton dimension of 'time' if (progvar(ivar)%numdims == 1) then if ((debug > 8) .and. do_output()) then write(*,*) write(*,*)'variable ',trim(progvar(ivar)%varname) write(*,*)'dimension 1 (i) ',progvar(ivar)%dimnames(1),progvar(ivar)%dimlens(1) endif SELECT CASE ( trim(progvar(ivar)%dimnames(1)) ) CASE ("gridcell") do i = 1, progvar(ivar)%dimlens(1) xi = grid1d_ixy(i) xj = grid1d_jxy(i) ! always unity if unstructured if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) endif indx = indx + 1 enddo CASE ("landunit") do i = 1, progvar(ivar)%dimlens(1) xi = land1d_ixy(i) xj = land1d_jxy(i) ! always unity if unstructured if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * land1d_wtxy(i) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * land1d_wtxy(i) endif indx = indx + 1 enddo CASE ("column") do i = 1, progvar(ivar)%dimlens(1) xi = cols1d_ixy(i) xj = cols1d_jxy(i) ! always unity if unstructured if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * cols1d_wtxy(i) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * cols1d_wtxy(i) endif indx = indx + 1 enddo CASE ("pft") do i = 1, progvar(ivar)%dimlens(1) xi = pfts1d_ixy(i) xj = pfts1d_jxy(i) ! always unity if unstructured if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * pfts1d_wtxy(i) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * pfts1d_wtxy(i) endif indx = indx + 1 enddo CASE DEFAULT write(string1,*)'(1d) unknown Dimension name '//trim(progvar(ivar)%dimnames(1))// & ' while trying to create metadata for '//trim(progvar(ivar)%varname) call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate) END SELECT elseif (progvar(ivar)%numdims == 2) then ! restart file variables may have 2 dimensions, one of which is layers. ! vector_history files may have 2 dimensions, one of which is time ! history file variables always have 3 dimension [time,lat,lon] ! In the ncdump output, dimension 2 is the leftmost dimension. ! Only dimension 2 matters for the weights. if ((debug > 8) .and. do_output()) then write(*,*) write(*,*)'variable ',trim(progvar(ivar)%varname) write(*,*)'dimension 1 (i) ',progvar(ivar)%dimnames(1),progvar(ivar)%dimlens(1) write(*,*)'dimension 2 (j) ',progvar(ivar)%dimnames(2),progvar(ivar)%dimlens(2) endif SELECT CASE ( trim(progvar(ivar)%dimnames(2)) ) CASE ("gridcell") if ((debug > 8) .and. do_output()) write(*,*)'length grid1d_ixy ',size(grid1d_ixy) do j = 1, progvar(ivar)%dimlens(2) xi = grid1d_ixy(j) xj = grid1d_jxy(j) ! always unity if unstructured do i = 1, progvar(ivar)%dimlens(1) if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) endif indx = indx + 1 enddo enddo CASE ("landunit") if ((debug > 8) .and. do_output()) write(*,*)'length land1d_ixy ',size(land1d_ixy) do j = 1, progvar(ivar)%dimlens(2) xi = land1d_ixy(j) xj = land1d_jxy(j) ! always unity if unstructured do i = 1, progvar(ivar)%dimlens(1) if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * land1d_wtxy(j) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * land1d_wtxy(j) endif indx = indx + 1 enddo enddo CASE ("column") ! Column is the only coordinate that has vertical levels. ! The vertical levels are fully defined by the levgrnd and ! levsno variables. levgrnd is static, levsno varies by column. progvar(ivar)%maxlevels = progvar(ivar)%dimlens(2) if ((debug > 8) .and. do_output()) write(*,*)'length cols1d_ixy ',size(cols1d_ixy) if ((debug > 8) .and. do_output()) write(*,*)'size zsno ',size(zsno,1), size(zsno,2) if ((debug > 8) .and. do_output()) write(*,*)'nlevsno ',nlevsno LANDCOLUMN : do j = 1, progvar(ivar)%dimlens(2) call fill_levels(progvar(ivar)%dimnames(1),j,progvar(ivar)%dimlens(1),levtot) xi = cols1d_ixy(j) xj = cols1d_jxy(j) ! always unity if unstructured VERTICAL : do i = 1, progvar(ivar)%dimlens(1) levels( indx) = levtot(i) if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * cols1d_wtxy(j) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * cols1d_wtxy(j) endif indx = indx + 1 enddo VERTICAL enddo LANDCOLUMN CASE ("pft") if ((debug > 8) .and. do_output()) write(*,*)'length pfts1d_ixy ',size(pfts1d_ixy) do j = 1, progvar(ivar)%dimlens(2) xi = pfts1d_ixy(j) xj = pfts1d_jxy(j) ! always unity if unstructured do i = 1, progvar(ivar)%dimlens(1) if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * pfts1d_wtxy(j) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * pfts1d_wtxy(j) endif indx = indx + 1 enddo enddo CASE ("time") ! The vector history files can have things 'pft,time' or 'column,time' SELECT CASE ( trim(progvar(ivar)%dimnames(1)) ) CASE ("pft") do i = 1, progvar(ivar)%dimlens(1) xi = pfts1d_ixy(i) xj = pfts1d_jxy(i) ! always unity if unstructured if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * pfts1d_wtxy(i) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * pfts1d_wtxy(i) endif indx = indx + 1 enddo CASE ("column") do i = 1, progvar(ivar)%dimlens(1) xi = cols1d_ixy(i) xj = cols1d_jxy(i) ! always unity if unstructured if (unstructured) then lonixy( indx) = xi latjxy( indx) = xi landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * cols1d_wtxy(i) else lonixy( indx) = xi latjxy( indx) = xj landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * cols1d_wtxy(i) endif indx = indx + 1 enddo CASE DEFAULT write(string1,*)'(2d) unknown first Dimension name '//trim(progvar(ivar)%dimnames(1))// & ' while trying to create metadata for '//trim(progvar(ivar)%varname) call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate) END SELECT CASE DEFAULT write(string1,*)'(2d) unknown Dimension name '//trim(progvar(ivar)%dimnames(2))// & ' while trying to create metadata for '//trim(progvar(ivar)%varname) call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate) END SELECT elseif (progvar(ivar)%numdims == 3) then ! restart file variables never have 3 dimensions ! vector_history variables may have 3 dimensions [time, lat, lon] ! history file variables always have 3 dimensions [time, lat, lon] ! exception is float H2OSOI(time, levgrnd, lat, lon) ... but we ! have access to restart file h2osoi_[liq,ice] if ((debug > 8) .and. do_output()) then write(*,*) write(*,*)'variable ',trim(progvar(ivar)%varname) write(*,*)'dimension 1 (i) ',progvar(ivar)%dimnames(1),progvar(ivar)%dimlens(1) write(*,*)'dimension 2 (j) ',progvar(ivar)%dimnames(2),progvar(ivar)%dimlens(2) write(*,*)'dimension 3 (k) ',progvar(ivar)%dimnames(3),progvar(ivar)%dimlens(3) endif ! Remember the order is reversed from ncdump to fortran ! The messages are in ncdump-order, checking is fortran-order if ((trim(progvar(ivar)%dimnames(1)) .ne. 'lon') .or. & (trim(progvar(ivar)%dimnames(2)) .ne. 'lat' ) .or. & (trim(progvar(ivar)%dimnames(3)) .ne. 'time' )) then write(string1,*)'3D variables must be [time,lat,lon] (as reported by ncdump)' write(string2,*)trim(progvar(ivar)%varname),' is ', & trim(progvar(ivar)%dimnames(3)), ' ', & trim(progvar(ivar)%dimnames(2)), ' ', & trim(progvar(ivar)%dimnames(1)) call error_handler(E_ERR,'static_init_model', string1, & source, revision, revdate, text2=string2) endif ! The get_var_3d() routine ensures there is only 1 timestep. ! So there is no need to loop over dimlens(3) do j = 1, progvar(ivar)%dimlens(2) ! time-invariant do i = 1, progvar(ivar)%dimlens(1) ! time-invariant lonixy( indx) = i latjxy( indx) = j landarea(indx) = AREA2D(i,j) * LANDFRAC2D(i,j) indx = indx + 1 enddo enddo else ! Cannot support 4D variables yet. ! These only occurr in 2D history files for variables with levels (and time) ! i.e. H2OSOI(time, levgrnd, lat, lon) ! ! You can get the same information from the restart file, usually. write(string1,*)'variables of rank ',progvar(ivar)%numdims,' are unsupported.' write(string2,*)trim(progvar(ivar)%varname),' is dimensioned ',& progvar(ivar)%dimlens(1:progvar(ivar)%numdims) call error_handler(E_ERR,'static_init_model', string1, source, revision, & revdate, text2=string2) endif indx = indx - 1 if (indx /= progvar(ivar)%indexN ) then write(string1,*)'variable ',trim(progvar(ivar)%varname), & ' is supposed to end at index ',progvar(ivar)%indexN write(string2,*)'it ends at index ',indx call error_handler(E_ERR,'static_init_model', string1, source, revision, & revdate, text2=string2) endif enddo end subroutine static_init_model subroutine end_model() !------------------------------------------------------------------ ! ! Does any shutdown and clean-up needed for model. Can be a NULL ! INTERFACE if the model has no need to clean up storage, etc. ! if ( .not. module_initialized ) call static_init_model if (unstructured) then deallocate(AREA1D, LANDFRAC1D) else deallocate(AREA2D, LANDFRAC2D) endif deallocate(LAT, LON, LEVGRND) deallocate(grid1d_ixy, grid1d_jxy) deallocate(land1d_ixy, land1d_jxy, land1d_wtxy) deallocate(cols1d_ixy, cols1d_jxy, cols1d_wtxy, cols1d_ityplun) deallocate(pfts1d_ixy, pfts1d_jxy, pfts1d_wtxy) deallocate(ens_mean) deallocate(lonixy, latjxy, landarea) end subroutine end_model subroutine init_time(time) !------------------------------------------------------------------ ! ! Companion interface to init_conditions. Returns a time that is somehow ! appropriate for starting up a long integration of the model. ! At present, this is only used if the namelist parameter ! start_from_restart is set to .false. in the program perfect_model_obs. ! If this option is not to be used in perfect_model_obs, or if no ! synthetic data experiments using perfect_model_obs are planned, ! this can be a NULL INTERFACE. type(time_type), intent(out) :: time if ( .not. module_initialized ) call static_init_model ! for now, just set to 0 time = set_time(0,0) end subroutine init_time subroutine init_conditions(x) !------------------------------------------------------------------ ! ! Returns a model state vector, x, that is some sort of appropriate ! initial condition for starting up a long integration of the model. ! At present, this is only used if the namelist parameter ! start_from_restart is set to .false. in the program perfect_model_obs. ! If this option is not to be used in perfect_model_obs, or if no ! synthetic data experiments using perfect_model_obs are planned, ! this can be a NULL INTERFACE. real(r8), intent(out) :: x(:) if ( .not. module_initialized ) call static_init_model x = 0.0_r8 end subroutine init_conditions function nc_write_model_atts( ncFileID ) result (ierr) !------------------------------------------------------------------ ! TJH -- Writes the model-specific attributes to a netCDF file. ! This includes coordinate variables and some metadata, but NOT ! the model state vector. ! ! assim_model_mod:init_diag_output uses information from the location_mod ! to define the location dimension and variable ID. All we need to do ! is query, verify, and fill ... ! ! Typical sequence for adding new dimensions,variables,attributes: ! NF90_OPEN ! open existing netCDF dataset ! NF90_redef ! put into define mode ! NF90_def_dim ! define additional dimensions (if any) ! NF90_def_var ! define variables: from name, type, and dims ! NF90_put_att ! assign attribute values ! NF90_ENDDEF ! end definitions: leave define mode ! NF90_put_var ! provide values for variable ! NF90_CLOSE ! close: save updated netCDF dataset integer, intent(in) :: ncFileID ! netCDF file identifier integer :: ierr ! return value of function integer :: nDimensions, nVariables, nAttributes, unlimitedDimID !---------------------------------------------------------------------- ! variables if we just blast out one long state vector !---------------------------------------------------------------------- integer :: StateVarDimID ! netCDF pointer to state variable dimension (model size) integer :: MemberDimID ! netCDF pointer to dimension of ensemble (ens_size) integer :: TimeDimID ! netCDF pointer to time dimension (unlimited) integer :: StateVarVarID ! netCDF pointer to state variable coordinate array integer :: StateVarID ! netCDF pointer to 3D [state,copy,time] array !---------------------------------------------------------------------- ! variables if we parse the state vector into prognostic variables. !---------------------------------------------------------------------- ! for the dimensions and coordinate variables integer :: nlonDimID integer :: nlatDimID integer :: gridcellDimID integer :: landunitDimID integer :: columnDimID integer :: pftDimID integer :: levgrndDimID integer :: levlakDimID integer :: levsnoDimID integer :: levsno1DimID integer :: levtotDimID integer :: numradDimID ! for the prognostic variables integer :: ivar, VarID !---------------------------------------------------------------------- ! local variables !---------------------------------------------------------------------- ! we are going to need these to record the creation date in the netCDF file. ! This is entirely optional, but nice. character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic character(len=5) :: crzone ! needed by F90 DATE_AND_TIME intrinsic integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic character(len=NF90_MAX_NAME) :: str1 character(len=NF90_MAX_NAME) :: varname integer, dimension(NF90_MAX_VAR_DIMS) :: mydimids integer :: i, myndims character(len=128) :: filename if ( .not. module_initialized ) call static_init_model ierr = -1 ! assume things go poorly !-------------------------------------------------------------------- ! we only have a netcdf handle here so we do not know the filename ! or the fortran unit number. but construct a string with at least ! the netcdf handle, so in case of error we can trace back to see ! which netcdf file is involved. !-------------------------------------------------------------------- write(filename,*) 'ncFileID', ncFileID !------------------------------------------------------------------------------- ! make sure ncFileID refers to an open netCDF file, ! and then put into define mode. !------------------------------------------------------------------------------- call nc_check(nf90_Inquire(ncFileID,nDimensions,nVariables,nAttributes,unlimitedDimID),& 'nc_write_model_atts', 'inquire '//trim(filename)) call nc_check(nf90_Redef(ncFileID),'nc_write_model_atts', 'redef '//trim(filename)) !------------------------------------------------------------------------------- ! We need the dimension ID for the number of copies/ensemble members, and ! we might as well check to make sure that Time is the Unlimited dimension. ! Our job is create the 'model size' dimension. !------------------------------------------------------------------------------- call nc_check(nf90_inq_dimid(ncid=ncFileID, name='copy', dimid=MemberDimID), & 'nc_write_model_atts', 'copy dimid '//trim(filename)) call nc_check(nf90_inq_dimid(ncid=ncFileID, name='time', dimid= TimeDimID), & 'nc_write_model_atts', 'time dimid '//trim(filename)) if ( TimeDimID /= unlimitedDimId ) then write(string1,*)'Time Dimension ID ',TimeDimID, & ' should equal Unlimited Dimension ID',unlimitedDimID call error_handler(E_ERR,'nc_write_model_atts', string1, source, revision, revdate) endif !------------------------------------------------------------------------------- ! Define the model size / state variable dimension / whatever ... !------------------------------------------------------------------------------- call nc_check(nf90_def_dim(ncid=ncFileID, name='StateVariable', len=model_size, & dimid = StateVarDimID),'nc_write_model_atts', 'state def_dim '//trim(filename)) !------------------------------------------------------------------------------- ! Write Global Attributes !------------------------------------------------------------------------------- call DATE_AND_TIME(crdate,crtime,crzone,values) write(str1,'(''YYYY MM DD HH MM SS = '',i4,5(1x,i2.2))') & values(1), values(2), values(3), values(5), values(6), values(7) call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'creation_date' ,str1 ), & 'nc_write_model_atts', 'creation put '//trim(filename)) call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model_source' ,source ), & 'nc_write_model_atts', 'source put '//trim(filename)) call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model_revision',revision), & 'nc_write_model_atts', 'revision put '//trim(filename)) call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model_revdate' ,revdate ), & 'nc_write_model_atts', 'revdate put '//trim(filename)) call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model', 'clm' ), & 'nc_write_model_atts', 'model put '//trim(filename)) !------------------------------------------------------------------------------- ! Here is the extensible part. The simplest scenario is to output the state vector, ! parsing the state vector into model-specific parts is complicated, and you need ! to know the geometry, the output variables (PS,U,V,T,Q,...) etc. We're skipping ! complicated part. !------------------------------------------------------------------------------- if ( output_state_vector ) then !---------------------------------------------------------------------------- ! Create a variable for the state vector !---------------------------------------------------------------------------- ! Define the state vector coordinate variable and some attributes. call nc_check(nf90_def_var(ncid=ncFileID,name='StateVariable', xtype=nf90_int, & dimids=StateVarDimID, varid=StateVarVarID), 'nc_write_model_atts', & 'statevariable def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID,StateVarVarID,'long_name','State Variable ID'),& 'nc_write_model_atts','statevariable long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, StateVarVarID, 'units','indexical'), & 'nc_write_model_atts', 'statevariable units '//trim(filename)) call nc_check(nf90_put_att(ncFileID,StateVarVarID,'valid_range',(/ 1,model_size /)),& 'nc_write_model_atts', 'statevariable valid_range '//trim(filename)) ! Define the actual (3D) state vector, which gets filled as time goes on ... call nc_check(nf90_def_var(ncid=ncFileID, name='state', xtype=nf90_real, & dimids=(/StateVarDimID,MemberDimID,unlimitedDimID/),varid=StateVarID),& 'nc_write_model_atts','state def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID,StateVarID,'long_name','model state or fcopy'),& 'nc_write_model_atts', 'state long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID,StateVarID,'_FillValue',MISSING_R8),& 'nc_write_model_atts', 'state FillValue '//trim(filename)) call nc_check(nf90_put_att(ncFileID,StateVarID,'missing_value',MISSING_R8),& 'nc_write_model_atts', 'state missing_value '//trim(filename)) ! Leave define mode so we can fill the coordinate variable. call nc_check(nf90_enddef(ncfileID),'nc_write_model_atts','state enddef '//trim(filename)) ! Fill the state variable coordinate variable call nc_check(nf90_put_var(ncFileID, StateVarVarID, (/ (i,i=1,model_size) /) ), & 'nc_write_model_atts', 'state put_var '//trim(filename)) else !---------------------------------------------------------------------------- ! We need to output the prognostic variables. !---------------------------------------------------------------------------- ! Define the new dimensions IDs !---------------------------------------------------------------------------- call nc_check(nf90_def_dim(ncid=ncFileID, name='lon', len = nlon, & dimid = nlonDimID),'nc_write_model_atts', 'lon def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='lat', len = nlat, & dimid = NlatDimID),'nc_write_model_atts', 'lat def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='gridcell', len = ngridcell, & dimid = gridcellDimID),'nc_write_model_atts', 'gridcell def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='landunit', len = nlandunit, & dimid = landunitDimID),'nc_write_model_atts', 'landunit def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='column', len = ncolumn, & dimid = columnDimID),'nc_write_model_atts', 'column def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='pft', len = npft, & dimid = pftDimID),'nc_write_model_atts', 'pft def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='levgrnd', len = nlevgrnd, & dimid = levgrndDimID),'nc_write_model_atts', 'levgrnd def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='levlak', len = nlevlak, & dimid = levlakDimID),'nc_write_model_atts', 'levlak def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='levsno', len = nlevsno, & dimid = levsnoDimID),'nc_write_model_atts', 'levsno def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='levsno1', len = nlevsno1, & dimid = levsno1DimID),'nc_write_model_atts', 'levsno1 def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='levtot', len = nlevtot, & dimid = levtotDimID),'nc_write_model_atts', 'levtot def_dim '//trim(filename)) call nc_check(nf90_def_dim(ncid=ncFileID, name='numrad', len = nnumrad, & dimid = numradDimID),'nc_write_model_atts', 'numrad def_dim '//trim(filename)) !---------------------------------------------------------------------------- ! Create the (empty) Coordinate Variables and the Attributes !---------------------------------------------------------------------------- !-------------------------------------------------------------------------------kyh03092016 !call nc_check(nf90_def_var(ncFileID,name='tran_can_v_grd', xtype=nf90_real, & ! dimids=(/ NlonDimID,nlatDimID,MemberDimID,TimeDimID /), varid=VarID),& ! 'nc_write_model_atts', 'tran_can_v_grd def_var '//trim(filename)) call nc_check(nf90_def_var(ncFileID,name='tran_can_v_grd', xtype=nf90_real, & dimids=(/ NlonDimID,nlatDimID,MemberDimID/), varid=VarID),& 'nc_write_model_atts', 'tran_can_v_grd def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'Vegetation transmissivity'), & 'nc_write_model_atts', 'tran_can_v_grd long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'units', '-'), & 'nc_write_model_atts', 'tran_can_v_grd units '//trim(filename)) !-------------------------------------------------------------------------------kyh03092016 ! Grid Longitudes call nc_check(nf90_def_var(ncFileID,name='lon', xtype=nf90_real, & dimids=(/ NlonDimID /), varid=VarID),& 'nc_write_model_atts', 'lon def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'coordinate longitude'), & 'nc_write_model_atts', 'lon long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'cartesian_axis', 'X'), & 'nc_write_model_atts', 'lon cartesian_axis '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'units', 'degrees_east'), & 'nc_write_model_atts', 'lon units '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'valid_range', (/ 0.0_r8, 360.0_r8 /)), & 'nc_write_model_atts', 'lon valid_range '//trim(filename)) ! Grid Latitudes call nc_check(nf90_def_var(ncFileID,name='lat', xtype=nf90_real, & dimids=(/ nlatDimID /), varid=VarID),& 'nc_write_model_atts', 'lat def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'coordinate latitude'), & 'nc_write_model_atts', 'lat long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'cartesian_axis', 'Y'), & 'nc_write_model_atts', 'lat cartesian_axis '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'units', 'degrees_north'), & 'nc_write_model_atts', 'lat units '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID,'valid_range',(/ -90.0_r8, 90.0_r8 /)), & 'nc_write_model_atts', 'lat valid_range '//trim(filename)) ! subsurface levels call nc_check(nf90_def_var(ncFileID,name='levgrnd', xtype=nf90_real, & dimids=(/ levgrndDimID /), varid=VarID),& 'nc_write_model_atts', 'levgrnd def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'coordinate soil levels'), & 'nc_write_model_atts', 'levgrnd long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'cartesian_axis', 'Z'), & 'nc_write_model_atts', 'levgrnd cartesian_axis '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'units', 'm'), & 'nc_write_model_atts', 'levgrnd units '//trim(filename)) ! grid cell areas if (unstructured) then call nc_check(nf90_def_var(ncFileID,name='area', xtype=nf90_real, & dimids=(/ NlonDimID /), varid=VarID),& 'nc_write_model_atts', 'area def_var '//trim(filename)) else call nc_check(nf90_def_var(ncFileID,name='area', xtype=nf90_real, & dimids=(/ NlonDimID,nlatDimID /), varid=VarID),& 'nc_write_model_atts', 'area def_var '//trim(filename)) endif call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'grid cell areas'), & 'nc_write_model_atts', 'area long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'units', 'km^2'), & 'nc_write_model_atts', 'area units '//trim(filename)) ! grid cell land fractions if (unstructured) then call nc_check(nf90_def_var(ncFileID,name='landfrac', xtype=nf90_real, & dimids=(/ NlonDimID /), varid=VarID),& 'nc_write_model_atts', 'landfrac def_var '//trim(filename)) else call nc_check(nf90_def_var(ncFileID,name='landfrac', xtype=nf90_real, & dimids=(/ NlonDimID,nlatDimID /), varid=VarID),& 'nc_write_model_atts', 'landfrac def_var '//trim(filename)) endif call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'land fraction'), & 'nc_write_model_atts', 'landfrac long_name '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'units', 'km^2'), & 'nc_write_model_atts', 'landfrac units '//trim(filename)) ! longitude grid index for each column call nc_check(nf90_def_var(ncFileID,name='cols1d_ixy', xtype=nf90_int, & dimids=(/ columnDimID /), varid=VarID),& 'nc_write_model_atts', 'cols1d_ixy def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & '2d longitude index of corresponding column'), & 'nc_write_model_atts', 'cols1d_ixy long_name '//trim(filename)) ! latitude grid index for each column call nc_check(nf90_def_var(ncFileID,name='cols1d_jxy', xtype=nf90_int, & dimids=(/ columnDimID /), varid=VarID),& 'nc_write_model_atts', 'cols1d_jxy def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & '2d latitude index of corresponding column'), & 'nc_write_model_atts', 'cols1d_jxy long_name '//trim(filename)) ! column weight relative to corresponding gridcell call nc_check(nf90_def_var(ncFileID,name='cols1d_wtxy', xtype=nf90_double, & dimids=(/ columnDimID /), varid=VarID),& 'nc_write_model_atts', 'cols1d_wtxy def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & 'column weight relative to corresponding gridcell'), & 'nc_write_model_atts', 'cols1d_wtxy long_name '//trim(filename)) ! land type to corresponding gridcell call nc_check(nf90_def_var(ncFileID,name='cols1d_ityplun', xtype=nf90_int, & dimids=(/ columnDimID /), varid=VarID),& 'nc_write_model_atts', 'cols1d_ityplun def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & 'column landunit type (vegetated,urban,lake,wetland or glacier)'), & 'nc_write_model_atts', 'cols1d_ityplun long_name '//trim(filename)) ! longitude grid index for each pft call nc_check(nf90_def_var(ncFileID,name='pfts1d_ixy', xtype=nf90_int, & dimids=(/ pftDimID /), varid=VarID),& 'nc_write_model_atts', 'pfts1d_ixy def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & '2d longitude index of corresponding column'), & 'nc_write_model_atts', 'pfts1d_ixy long_name '//trim(filename)) ! latitude grid index for each pft call nc_check(nf90_def_var(ncFileID,name='pfts1d_jxy', xtype=nf90_int, & dimids=(/ pftDimID /), varid=VarID),& 'nc_write_model_atts', 'pfts1d_jxy def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & '2d latitude index of corresponding column'), & 'nc_write_model_atts', 'pfts1d_jxy long_name '//trim(filename)) ! pft weight relative to corresponding gridcell call nc_check(nf90_def_var(ncFileID,name='pfts1d_wtxy', xtype=nf90_double, & dimids=(/ pftDimID /), varid=VarID),& 'nc_write_model_atts', 'pfts1d_wtxy def_var '//trim(filename)) call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', & 'pft weight relative to corresponding gridcell'), & 'nc_write_model_atts', 'pfts1d_wtxy long_name '//trim(filename)) !---------------------------------------------------------------------------- ! Create the (empty) Prognostic Variables and the Attributes !---------------------------------------------------------------------------- do ivar=1, nfields varname = trim(progvar(ivar)%varname) string1 = trim(filename)//' '//trim(varname) ! match shape of the variable to the dimension IDs call define_var_dims(ivar, ncFileID, MemberDimID, unlimitedDimID, myndims, mydimids) ! define the variable and set the attributes call nc_check(nf90_def_var(ncid=ncFileID, name=trim(varname), & xtype=progvar(ivar)%xtype, & dimids = mydimids(1:myndims), varid=VarID),& 'nc_write_model_atts', trim(string1)//' def_var' ) call nc_check(nf90_put_att(ncFileID, VarID, & 'long_name', trim(progvar(ivar)%long_name)), & 'nc_write_model_atts', trim(string1)//' put_att long_name' ) call nc_check(nf90_put_att(ncFileID, VarID, & 'DART_kind', trim(progvar(ivar)%kind_string)), & 'nc_write_model_atts', trim(string1)//' put_att dart_kind' ) call nc_check(nf90_put_att(ncFileID, VarID, & 'units', trim(progvar(ivar)%units)), & 'nc_write_model_atts', trim(string1)//' put_att units' ) ! Preserve the original missing_value/_FillValue code. if ( progvar(ivar)%xtype == NF90_INT ) then call nc_check(nf90_put_att(ncFileID, VarID, & 'missing_value', progvar(ivar)%spvalINT), & 'nc_write_model_atts', trim(string1)//' put_att missing_value' ) call nc_check(nf90_put_att(ncFileID, VarID, & '_FillValue', progvar(ivar)%spvalINT), & 'nc_write_model_atts', trim(string1)//' put_att _FillValue' ) elseif ( progvar(ivar)%xtype == NF90_FLOAT ) then call nc_check(nf90_put_att(ncFileID, VarID, & 'missing_value', progvar(ivar)%spvalR4), & 'nc_write_model_atts', trim(string1)//' put_att missing_value' ) call nc_check(nf90_put_att(ncFileID, VarID, & '_FillValue', progvar(ivar)%spvalR4), & 'nc_write_model_atts', trim(string1)//' put_att _FillValue' ) elseif ( progvar(ivar)%xtype == NF90_DOUBLE ) then call nc_check(nf90_put_att(ncFileID, VarID, & 'missing_value', progvar(ivar)%spvalR8), & 'nc_write_model_atts', trim(string1)//' put_att missing_value' ) call nc_check(nf90_put_att(ncFileID, VarID, & '_FillValue', progvar(ivar)%spvalR8), & 'nc_write_model_atts', trim(string1)//' put_att _FillValue' ) endif enddo !---------------------------------------------------------------------------- ! Finished with dimension/variable definitions, must end 'define' mode to fill. !---------------------------------------------------------------------------- call nc_check(nf90_enddef(ncfileID), 'prognostic enddef '//trim(filename)) !---------------------------------------------------------------------------- ! Fill the coordinate variables !---------------------------------------------------------------------------- call nc_check(nf90_inq_varid(ncFileID, 'lon', VarID), & 'nc_write_model_atts', 'put_var lon '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, lon ), & 'nc_write_model_atts', 'lon put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'lat', VarID), & 'nc_write_model_atts', 'put_var lat '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, lat ), & 'nc_write_model_atts', 'lat put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'levgrnd', VarID), & 'nc_write_model_atts', 'put_var levgrnd '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, LEVGRND ), & 'nc_write_model_atts', 'levgrnd put_var '//trim(filename)) ! AREA can be 1D or 2D call nc_check(nf90_inq_varid(ncFileID, 'area', VarID), & 'nc_write_model_atts', 'put_var area '//trim(filename)) if (unstructured) then call nc_check(nf90_put_var(ncFileID, VarID, AREA1D ), & 'nc_write_model_atts', 'area put_var '//trim(filename)) else call nc_check(nf90_put_var(ncFileID, VarID, AREA2D ), & 'nc_write_model_atts', 'area put_var '//trim(filename)) endif ! LANDFRAC can be 1D or 2D call nc_check(nf90_inq_varid(ncFileID, 'landfrac', VarID), & 'nc_write_model_atts', 'put_var landfrac '//trim(filename)) if (unstructured) then call nc_check(nf90_put_var(ncFileID, VarID, LANDFRAC1D ), & 'nc_write_model_atts', 'landfrac put_var '//trim(filename)) else call nc_check(nf90_put_var(ncFileID, VarID, LANDFRAC2D ), & 'nc_write_model_atts', 'landfrac put_var '//trim(filename)) endif call nc_check(nf90_inq_varid(ncFileID, 'cols1d_ixy', VarID), & 'nc_write_model_atts', 'put_var cols1d_ixy '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, cols1d_ixy ), & 'nc_write_model_atts', 'cols1d_ixy put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'cols1d_jxy', VarID), & 'nc_write_model_atts', 'put_var cols1d_jxy '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, cols1d_jxy ), & 'nc_write_model_atts', 'cols1d_jxy put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'cols1d_wtxy', VarID), & 'nc_write_model_atts', 'put_var cols1d_wtxy '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, cols1d_wtxy ), & 'nc_write_model_atts', 'cols1d_wtxy put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'cols1d_ityplun', VarID), & 'nc_write_model_atts', 'put_var cols1d_ityplun '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, cols1d_ityplun ), & 'nc_write_model_atts', 'cols1d_ityplun put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'pfts1d_ixy', VarID), & 'nc_write_model_atts', 'put_var pfts1d_ixy '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, pfts1d_ixy ), & 'nc_write_model_atts', 'pfts1d_ixy put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'pfts1d_jxy', VarID), & 'nc_write_model_atts', 'put_var pfts1d_jxy '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, pfts1d_jxy ), & 'nc_write_model_atts', 'pfts1d_jxy put_var '//trim(filename)) call nc_check(nf90_inq_varid(ncFileID, 'pfts1d_wtxy', VarID), & 'nc_write_model_atts', 'put_var pfts1d_wtxy '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, pfts1d_wtxy ), & 'nc_write_model_atts', 'pfts1d_wtxy put_var '//trim(filename)) endif !------------------------------------------------------------------------------- ! Fill the variables we can !------------------------------------------------------------------------------- ! none !------------------------------------------------------------------------------- ! Flush the buffer and leave netCDF file open !------------------------------------------------------------------------------- call nc_check(nf90_sync(ncFileID), 'nc_write_model_atts', 'atts sync') ierr = 0 ! If we got here, things went well. end function nc_write_model_atts function nc_write_model_vars( ncFileID, state_vec, copyindex, timeindex ) result (ierr) !------------------------------------------------------------------ ! Writes the model variables to a netCDF file. ! ! All errors are fatal, so the ! return code is always '0 == normal', since the fatal errors stop execution. ! ! assim_model_mod:init_diag_output uses information from the location_mod ! to define the location dimension and variable ID. All we need to do ! is query, verify, and fill ... ! ! Typical sequence for adding new dimensions,variables,attributes: ! NF90_OPEN ! open existing netCDF dataset ! NF90_redef ! put into define mode ! NF90_def_dim ! define additional dimensions (if any) ! NF90_def_var ! define variables: from name, type, and dims ! NF90_put_att ! assign attribute values ! NF90_ENDDEF ! end definitions: leave define mode ! NF90_put_var ! provide values for variable ! NF90_CLOSE ! close: save updated netCDF dataset integer, intent(in) :: ncFileID ! netCDF file identifier real(r8), dimension(:), intent(in) :: state_vec integer, intent(in) :: copyindex integer, intent(in) :: timeindex integer :: ierr ! return value of function integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, ncstart, nccount character(len=NF90_MAX_NAME) :: varname integer :: i, ivar, VarID, ncNdims, dimlen integer :: TimeDimID, CopyDimID real(r8), allocatable, dimension(:) :: data_1d_array real(r8), allocatable, dimension(:,:) :: data_2d_array real(r8), allocatable, dimension(:) :: data_1d_array_prior !kyh11142014 real(r8), allocatable, dimension(:,:) :: data_2d_array_prior !kyh11142014 character(len=128) :: filename if ( .not. module_initialized ) call static_init_model ierr = -1 ! assume things go poorly !-------------------------------------------------------------------- ! we only have a netcdf handle here so we do not know the filename ! or the fortran unit number. but construct a string with at least ! the netcdf handle, so in case of error we can trace back to see ! which netcdf file is involved. !-------------------------------------------------------------------- write(filename,*) 'ncFileID', ncFileID !------------------------------------------------------------------------------- ! make sure ncFileID refers to an open netCDF file, !------------------------------------------------------------------------------- call nc_check(nf90_inq_dimid(ncFileID, 'copy', dimid=CopyDimID), & 'nc_write_model_vars', 'inq_dimid copy '//trim(filename)) call nc_check(nf90_inq_dimid(ncFileID, 'time', dimid=TimeDimID), & 'nc_write_model_vars', 'inq_dimid time '//trim(filename)) if ( output_state_vector ) then call nc_check(NF90_inq_varid(ncFileID, 'state', VarID), & 'nc_write_model_vars', 'state inq_varid '//trim(filename)) call nc_check(NF90_put_var(ncFileID,VarID,state_vec,start=(/1,copyindex,timeindex/)),& 'nc_write_model_vars', 'state put_var '//trim(filename)) else !--------------------------------------------------------------------------------------kyh03092016 call nc_check(NF90_inq_varid(ncFileID, 'tran_can_v_grd', VarID), & 'nc_write_model_vars', 'tran_can_v_grd inq_varid '//trim(filename)) !call nc_check(NF90_put_var(ncFileID,VarID,tran_can_v_grd,start=(/1,copyindex,timeindex/)),& ! 'nc_write_model_vars', 'tran_can_v_grd put_var '//trim(filename)) call nc_check(NF90_put_var(ncFileID,VarID,tran_can_v_grd),& 'nc_write_model_vars', 'tran_can_v_grd put_var '//trim(filename)) ncstart(1:3) = (/ 1, 1, copyindex /) nccount(1:3) = (/ 288, 192, 1 /) call nc_check(NF90_put_var(ncFileID, VarID, tran_can_v_grd(:,:,copyindex), & start = ncstart(1:3), count=nccount(1:3)), & 'nc_write_model_vars', 'tran_can_v_grd put_var '//trim(filename)) ncstart(1:4) = (/ 1, 1, copyindex, timeindex /) nccount(1:4) = (/ 288, 192, 1, 1 /) call nc_check(NF90_put_var(ncFileID, VarID, tran_can_v_grd(:,:,copyindex), & start = ncstart(1:4), count=nccount(1:4)), & 'nc_write_model_vars', 'tran_can_v_grd put_var '//trim(filename)) !--------------------------------------------------------------------------------------kyh03092016 !---------------------------------------------------------------------------- ! We need to process the prognostic variables. !---------------------------------------------------------------------------- do ivar = 1,nfields ! Very similar to loop in sv_to_restart_file varname = trim(progvar(ivar)%varname) string2 = trim(filename)//' '//trim(varname) ! Ensure netCDF variable is conformable with progvar quantity. ! The TIME and Copy dimensions are intentionally not queried ! by looping over the dimensions stored in the progvar type. call nc_check(nf90_inq_varid(ncFileID, varname, VarID), & 'nc_write_model_vars', 'inq_varid '//trim(string2)) call nc_check(nf90_inquire_variable(ncFileID,VarID,dimids=dimIDs,ndims=ncNdims), & 'nc_write_model_vars', 'inquire '//trim(string2)) ncstart = 1 ! These are arrays, actually nccount = 1 DimCheck : do i = 1,progvar(ivar)%numdims write(string1,'(a,i2,A)') 'inquire dimension ',i,trim(string2) call nc_check(nf90_inquire_dimension(ncFileID, dimIDs(i), len=dimlen), & 'nc_write_model_vars', trim(string1)) if (progvar(ivar)%dimnames(i) == 'time') cycle DimCheck if ( dimlen /= progvar(ivar)%dimlens(i) ) then write(string1,*)trim(string2),' dim/dimlen ',i,dimlen, & ' not ',progvar(ivar)%dimlens(i) write(string2,*)' but it should be.' call error_handler(E_ERR, 'nc_write_model_vars', trim(string1), & source, revision, revdate, text2=trim(string2)) endif nccount(i) = dimlen enddo DimCheck where(dimIDs == CopyDimID) ncstart = copyindex where(dimIDs == CopyDimID) nccount = 1 where(dimIDs == TimeDimID) ncstart = timeindex where(dimIDs == TimeDimID) nccount = 1 if ((debug > 0) .and. do_output()) then write(*,*)'nc_write_model_vars '//trim(varname)//' start is ',ncstart(1:ncNdims) write(*,*)'nc_write_model_vars '//trim(varname)//' count is ',nccount(1:ncNdims) endif ! Since 'time' is a singleton dimension, we can use the same logic ! as if it the variable had one less dimension. if ( (progvar(ivar)%numdims == 1) .or. & ((progvar(ivar)%numdims == 2) .and. & (progvar(ivar)%dimnames(2) == 'time')) )then if ( ncNdims /= 3 ) then write(string1,*)trim(varname),' no room for copy,time dimensions.' write(string2,*)'netcdf file should have 3 dimensions, has ',ncNdims call error_handler(E_ERR, 'nc_write_model_vars', string1, & source, revision, revdate, text2=string2) endif allocate(data_1d_array( progvar(ivar)%dimlens(1) ) ) allocate(data_1d_array_prior( progvar(ivar)%dimlens(1) ) ) !kyh11142014 !call vector_to_prog_var(state_vec, ivar, data_1d_array) !original call vector_to_prog_var(state_vec, ivar, data_1d_array, data_1d_array_prior) !kyh11142014 call nc_check(nf90_put_var(ncFileID, VarID, data_1d_array, & start = ncstart(1:ncNdims), count=nccount(1:ncNdims)), & 'nc_write_model_vars', 'put_var '//trim(string2)) deallocate(data_1d_array) deallocate(data_1d_array_prior) !kyh11142014 elseif ( (progvar(ivar)%numdims == 2) .or. & ((progvar(ivar)%numdims == 3) .and. & (progvar(ivar)%dimnames(3) == 'time')) )then if ( ncNdims /= 4 ) then write(string1,*)trim(varname),' no room for copy,time dimensions.' write(string2,*)'netcdf file should have 4 dimensions, has ',ncNdims call error_handler(E_ERR, 'nc_write_model_vars', string1, & source, revision, revdate, text2=string2) endif allocate(data_2d_array( progvar(ivar)%dimlens(1), & progvar(ivar)%dimlens(2) )) allocate(data_2d_array_prior( progvar(ivar)%dimlens(1), & progvar(ivar)%dimlens(2) )) !kyh11142014 !call vector_to_prog_var(state_vec, ivar, data_2d_array) !original call vector_to_prog_var(state_vec, ivar, data_2d_array, data_2d_array_prior) !kyh11142014 call nc_check(nf90_put_var(ncFileID, VarID, data_2d_array, & start = ncstart(1:ncNdims), count=nccount(1:ncNdims)), & 'nc_write_model_vars', 'put_var '//trim(string2)) deallocate(data_2d_array) deallocate(data_2d_array_prior) !kyh11142014 else write(string1,*)'do not know how to handle CLM variables with more than 2 dimensions' write(string2,*)trim(progvar(ivar)%varname),'has shape', & progvar(ivar)%dimlens(1:progvar(ivar)%numdims) call error_handler(E_ERR,'nc_write_model_vars',string1,source,revision,revdate) endif enddo endif !------------------------------------------------------------------------------- ! Flush the buffer and leave netCDF file open !------------------------------------------------------------------------------- call nc_check(nf90_sync(ncFileID), 'nc_write_model_vars', 'sync '//trim(filename)) ierr = 0 ! If we got here, things went well. end function nc_write_model_vars subroutine pert_model_state(state, pert_state, interf_provided) !------------------------------------------------------------------ ! ! Perturbs a single model state for generating initial ensembles. ! This (required interface) is unsupported in CLM and any attempt ! to use it will cause DART to terminate. ! ! So far, we have generated intial ensembles by taking a single ! instance and replicating it N times - and pairing each of the ! instances with a unique atmospheric forcing file and integrating ! for some period of time till the ensemble demonstrates sufficient ! spread. This is an area of active research. ! ! The naieve approach does not work -- it generates negative ! snow cover fractions, for example. Must check for out-of-range ! values specific to each type. ! The WRF model mod has something that might be useful. real(r8), intent(in) :: state(:) real(r8), intent(out) :: pert_state(:) logical, intent(out) :: interf_provided integer :: i logical, save :: random_seq_init = .false. if ( .not. module_initialized ) call static_init_model call error_handler(E_ERR,'pert_model_state', & 'CLM cannot be started from a single vector', & source, revision, revdate, & text2='see comments in clm/model_mod.f90::pert_model_state()') interf_provided = .true. ! Initialize my random number sequence if(.not. random_seq_init) then call init_random_seq(random_seq, my_task_id()) random_seq_init = .true. endif ! This does keep the compiler error messages down, but this ! section of code must never be reached. do i=1,size(state) pert_state(i) = random_gaussian(random_seq, state(i), & model_perturbation_amplitude) enddo end subroutine pert_model_state subroutine ens_mean_for_model(filter_ens_mean) !------------------------------------------------------------------ ! If needed by the model interface, this is the current mean ! for all state vector items across all ensembles. real(r8), intent(in) :: filter_ens_mean(:) if ( .not. module_initialized ) call static_init_model ens_mean = filter_ens_mean end subroutine ens_mean_for_model !================================================================== ! The remaining PUBLIC interfaces come next !================================================================== subroutine clm_to_dart_state_vector(state_vector, restart_time) !------------------------------------------------------------------ ! Reads the current time and state variables from a clm restart ! file and packs them into a dart state vector. This better happen ! in the same fashion as the metadata arrays are built. real(r8), intent(out) :: state_vector(:) type(time_type), intent(out) :: restart_time ! temp space to hold data while we are reading it integer :: i, j, ni, nj, ivar, indx, numsnowlevels integer, allocatable, dimension(:) :: snlsno real(r8), allocatable, dimension(:) :: data_1d_array real(r8), allocatable, dimension(:,:) :: data_2d_array real(r8), allocatable, dimension(:,:,:) :: data_3d_array integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs character(len=NF90_MAX_NAME) :: varname integer :: io, TimeDimID, VarID, ncNdims, dimlen integer :: ncid character(len=256) :: myerrorstring if ( .not. module_initialized ) call static_init_model state_vector = MISSING_R8 ! Must check anything with a dimension of 'levtot' or 'levsno' and manually ! set the values to DART missing. If only it were that easy ... ! ! The treatment of snow-related variables is complicated. ! The SNLSNO variable defines the number of snow layers with valid values. ! HOWEVER, if the snow depth is < 0.01 m, the snow is not represented by a layer, ! so the SNLSNO(i) is zero even though there is a trace of snow. ! Even a trace amount of snow results in some sort of snow cover fraction. ! ! Lakes are treated differently. ! The SNLSNO(i) is always zero, even though there is snow. ! The snow over lakes is wholly contained in the bulk formulation variables ! as opposed to the snow layer variables. ! Some things must be gotten explicitly from the restart file - ONCE. ! number of snow layers ! time of restart file allocate(snlsno(Ncolumn)) call nc_check(nf90_open(trim(clm_restart_filename), NF90_NOWRITE, ncid), & 'clm_to_dart_state_vector', 'open SNLSNO'//clm_restart_filename) call nc_check(nf90_inq_varid(ncid,'SNLSNO', VarID), & 'clm_to_dart_state_vector', 'inq_varid SNLSNO'//clm_restart_filename) call nc_check(nf90_get_var(ncid, VarID, snlsno), & 'clm_to_dart_state_vector', 'get_var SNLSNO'//clm_restart_filename) restart_time = get_state_time(ncid) if (do_output()) call print_time(restart_time,'time in restart file '//clm_restart_filename) if (do_output()) call print_date(restart_time,'date in restart file '//clm_restart_filename) call nc_check(nf90_close(ncid),'clm_to_dart_state_vector','close '//clm_restart_filename) ! Start counting and filling the state vector one item at a time, ! repacking the Nd arrays into a single 1d list of numbers. do ivar=1, nfields varname = trim(progvar(ivar)%varname) myerrorstring = trim(progvar(ivar)%origin)//' '//trim(progvar(ivar)%varname) call nc_check(nf90_open(trim(progvar(ivar)%origin), NF90_NOWRITE, ncid), & 'clm_to_dart_state_vector','open '//trim(myerrorstring)) ! File is not required to have a time dimension io = nf90_inq_dimid(ncid, 'time', TimeDimID) if (io /= NF90_NOERR) TimeDimID = MISSING_I call nc_check(nf90_inq_varid(ncid, varname, VarID), & 'clm_to_dart_state_vector', 'inq_varid '//trim(myerrorstring)) call nc_check(nf90_inquire_variable(ncid,VarID,dimids=dimIDs,ndims=ncNdims), & 'clm_to_dart_state_vector', 'inquire '//trim(myerrorstring)) ! Check the rank of the variable if ( ncNdims /= progvar(ivar)%numdims ) then write(string1, *) 'netCDF rank of '//trim(varname)//' does not match derived type knowledge' write(string2, *) 'netCDF rank is ',ncNdims,' expected ',progvar(ivar)%numdims call error_handler(E_ERR,'clm_to_dart_state_vector', string1, & source,revision,revdate,text2=string2) endif ! Check the shape of the variable do i = 1,progvar(ivar)%numdims write(string1,'(''inquire dimension'',i2,A)') i,trim(myerrorstring) call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), & 'clm_to_dart_state_vector', string1) ! Time dimension will be unity in progvar, but not necessarily ! in origin file. We only want a single matching time. ! static_init_model() only reserves space for a single time. if ( dimIDs(i) == TimeDimID ) dimlen = 1 if ( dimlen /= progvar(ivar)%dimlens(i) ) then write(string1,*) trim(myerrorstring),' dim/dimlen ',i,dimlen, & ' not ',progvar(ivar)%dimlens(i) call error_handler(E_ERR,'clm_to_dart_state_vector',string1,source,revision,revdate) endif enddo ! Pack the variable into the DART state vector ! Could/should fill metadata arrays at the same time ... ! As of 24 Aug 2011, CLM was not consistent about using a single fill_value ! or missing value, and the restart files didn't use the right attributes anyway ... ! (bugzilla report 1401) indx = progvar(ivar)%index1 if (ncNdims == 1) then ni = progvar(ivar)%dimlens(1) allocate(data_1d_array(ni)) call DART_get_var(ncid, varname, data_1d_array) do i = 1, ni state_vector(indx) = data_1d_array(i) indx = indx + 1 enddo deallocate(data_1d_array) elseif (ncNdims == 2) then ni = progvar(ivar)%dimlens(1) nj = progvar(ivar)%dimlens(2) allocate(data_2d_array(ni, nj)) call DART_get_var(ncid, varname, data_2d_array) ! README: The values in unused snow layers must be assumed to be ! indeterminate. If the layer is not in use, fill with a missing value. ! (levsno,column) and (levtot,column) variables may be treated identically. ! abs(snlsno(j)) defines the number of valid levels in each column - ! even over lakes. Lakes use a 'bulk' formula, so all the pertinent ! values are in the 1D variables, SNOWDP and frac_sno. ! FIXME: Question, what happens to unused levels below ground? Are those ! values 'special'? if ( (trim(progvar(ivar)%dimnames(1)) == 'levsno') .and. & (trim(progvar(ivar)%dimnames(2)) == 'column') ) then do j = 1, nj ! loop over columns numsnowlevels = abs(snlsno(j)) do i = 1, Nlevsno - numsnowlevels ! loop over layers data_2d_array(i,j) = MISSING_R8 enddo enddo elseif ( (trim(progvar(ivar)%dimnames(1)) == 'levtot') .and. & (trim(progvar(ivar)%dimnames(2)) == 'column') ) then do j = 1, nj ! loop over columns numsnowlevels = abs(snlsno(j)) do i = 1, Nlevsno - numsnowlevels ! loop over layers data_2d_array(i,j) = MISSING_R8 enddo enddo endif ! Block of checks that will hopefully be corrected in the ! core CLM code. There are some indeterminate values being ! used instead of the missing_value code - and even then, ! the missing_value code is not reliably implemented. if (progvar(ivar)%varname == 'T_SOISNO') then where(data_2d_array < 1.0_r8) data_2d_array = MISSING_R8 do j = 1,nj ! T_SOISNO has missing data in lake columns if (cols1d_ityplun(j) == LAKE) then ! write(*,*)'Found a lake column resetting the following:' ! write(*,*)data_2d_array(:,j) data_2d_array(:,j) = MISSING_R8 endif enddo endif if ((progvar(ivar)%varname == 'H2OSOI_LIQ') .or. & (progvar(ivar)%varname == 'H2OSOI_ICE')) then where(data_2d_array < 0.0_r8) data_2d_array = MISSING_R8 do j = 1,nj ! missing data in lake columns if (cols1d_ityplun(j) == LAKE) then data_2d_array(:,j) = MISSING_R8 endif enddo endif ! Finally, pack it into the DART state vector. do j = 1, nj do i = 1, ni state_vector(indx) = data_2d_array(i, j) indx = indx + 1 enddo enddo deallocate(data_2d_array) elseif (ncNdims == 3) then ! restart file variables never have 3 dimensions ! vector_history variables may have 3 dimensions [time, lat, lon] ! history file variables always have 3 dimensions [time, lat, lon] ! exception is float H2OSOI(time, levgrnd, lat, lon) ... but we ! have access to restart file h2osoi_[liq,ice] if ( (trim(progvar(ivar)%dimnames(1)) == 'lon') .and. & (trim(progvar(ivar)%dimnames(2)) == 'lat') .and. & (trim(progvar(ivar)%dimnames(3)) == 'time') ) then ni = progvar(ivar)%dimlens(1) nj = progvar(ivar)%dimlens(2) ! nk = progvar(ivar)%dimlens(3) not needed ... time is always a singleton allocate(data_3d_array(ni, nj, 1)) call DART_get_var(ncid, varname, data_3d_array) ! In the CLM history files, the _missing_value_ flag seems to be ! applied correctly for PBOT, TBOT ... so there is no need for the ! extra processing that is present in the previous loops. do j = 1, nj do i = 1, ni state_vector(indx) = data_3d_array(i, j, 1) indx = indx + 1 enddo enddo deallocate(data_3d_array) else write(string1, *) '3D variable unexpected shape -- only support nlon, nlat, time(=1)' write(string2, *) 'variable [',trim(progvar(ivar)%varname),']' write(string3, *) 'file [',trim(progvar(ivar)%origin),']' call error_handler(E_ERR,'clm_to_dart_state_vector', string1, & source, revision, revdate, text2=string2, text3=string3) endif else write(string1, *) 'no support for data array of dimension ', ncNdims write(string2, *) 'variable [',trim(progvar(ivar)%varname),']' write(string3, *) 'file [',trim(progvar(ivar)%origin),']' call error_handler(E_ERR,'clm_to_dart_state_vector', string1, & source, revision, revdate, text2=string2, text3=string3) endif indx = indx - 1 if ( indx /= progvar(ivar)%indexN ) then write(string1, *)'Variable '//trim(varname)//' filled wrong.' write(string2, *)'Should have ended at ',progvar(ivar)%indexN,' actually ended at ',indx call error_handler(E_ERR,'clm_to_dart_state_vector', string1, & source,revision,revdate,text2=string2) endif call nc_check(nf90_close(ncid),'clm_to_dart_state_vector','close '//progvar(ivar)%origin) ncid = 0 enddo deallocate(snlsno) end subroutine clm_to_dart_state_vector subroutine sv_to_restart_file(state_vector, filename, dart_time) !------------------------------------------------------------------ ! Writes the current time and state variables from a dart state ! vector (1d array) into a clm netcdf restart file. ! real(r8), intent(in) :: state_vector(:) character(len=*), intent(in) :: filename type(time_type), intent(in) :: dart_time ! temp space to hold data while we are writing it integer :: i, ni, nj, ivar integer :: j,c !kyh11142014 integer, allocatable, dimension(:) :: snlsno !kyh11142014 real(r8), allocatable, dimension(:) :: data_1d_array real(r8), allocatable, dimension(:,:) :: data_2d_array real(r8), allocatable, dimension(:) :: data_1d_array_prior !kyh11142014 real(r8), allocatable, dimension(:,:) :: data_2d_array_prior !kyh11142014 real(r8), allocatable, dimension(:) :: h2osno_pr,h2osno_po,snowdp_pr,snowdp_po, gain_snowdp !kyh11142014 real(r8), allocatable, dimension(:,:) :: gain_h2osno_l, gain_dzsno, gain_h2oliq_l,& gain_h2oice_l, dzsno_pr, dzsno_po, h2oliq_pr, h2oliq_po, h2oice_pr, & h2oice_po !kyh11142014 real(r8) :: wt_liq_l, snowden, wt_swe_l, reduce_factor !kyh11142014 real(r8), parameter :: NaN = O'0777610000000000000000' !kyh11142014 integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs character(len=NF90_MAX_NAME) :: varname integer :: VarID, ncNdims, dimlen integer :: ncFileID type(time_type) :: file_time if ( .not. module_initialized ) call static_init_model ! Check that the output file exists ... if ( .not. file_exist(filename) ) then write(string1,*) 'cannot open file ', trim(filename),' for writing.' call error_handler(E_ERR,'sv_to_restart_file',string1,source,revision,revdate) endif call nc_check(nf90_open(trim(filename), NF90_WRITE, ncFileID), & 'sv_to_restart_file','open '//trim(filename)) ! make sure the time in the file is the same as the time on the data ! we are trying to insert. we are only updating part of the contents ! of the clm restart file, and state vector contents from a different ! time won't be consistent with the rest of the file. file_time = get_state_time(ncFileID) if ( file_time /= dart_time ) then call print_time(dart_time,'DART current time',logfileunit) call print_time(file_time,'clm current time',logfileunit) call print_time(dart_time,'DART current time') call print_time(file_time,'clm current time') write(string1,*)trim(filename),' current time /= model time. FATAL error.' call error_handler(E_ERR,'sv_to_restart_file',string1,source,revision,revdate) endif if (do_output()) call print_time(file_time,'time of restart file '//trim(filename)) if (do_output()) call print_date(file_time,'date of restart file '//trim(filename)) ! The DART prognostic variables are only defined for a single time. ! We already checked the assumption that variables are xy2d or xyz3d ... ! IF the netCDF variable has a TIME dimension, it must be the last dimension, ! and we need to read the LAST timestep and effectively squeeze out the ! singleton dimension when we stuff it into the DART state vector. ! The snow water equivalent (H2OSNO) cannot be zero since H2OSNO is used to calculate the ! bulk snow density, which in turn is a parameter in the equation of Snow Cover Fraction. ! In order to avoid the negative values of H2OSNO produced by DART, I added some "if" conditions ! to set the value of H2OSNO back to the value before assimilation if negative value is found. !------------------------------------------------------------------------------------------------------kyh05072014 !!========================================================================= ! Test codes written by Fei. I intend to send DZSNO,H2OSOI_ICE, H2OSOI_LIQ, ! H2OSNO and SNOWDP back to CLM4. But only H2OSNO is used to be updated ! directly by the filter. The update of other variables will be calculated ! based on their physical relationships. ! SnowDensity(layer,column) = (H2OSOI_LIQ(layer,column)+H2OSOI_ICE(layer,column)) ! /DZSNO(layer,column) ! GainH2OSNO_l(layer,column) = GainH2OSNO(column)*wt_swe(layer,column) ! wt_swe(layer,column) = H2OSOI_LIQ(layer,column)+H2OSOI_ICE(layer,column)) ! /H2OSNO(column) ! GainDZSNO(layer,column) = GainH2OSNO_l(layer,column)/SnowDensity(layer,column) ! wt_liq(layer,column) = H2OSOI_LIQ(layer,column)/(H2OSOI_LIQ(layer,column)+ ! H2OSOI_ICE(layer,column)) ! wt_ice(layer,column) = 1-wt_liq(layer,column) ! GainH2O_LIQ(layer,column) = GainH2OSNO_l(layer,column)*wt_liq(layer,column) ! GainH2O_ICE(layer,column) = GainH2OSNO_l(layer,column)*(1-wt_liq(layer,column)) ! GainSNOWDP(column) = sum(GainDZSNO(layer,column)) ! SNOWDP_update(column) = SNOWDP(column)+GainSNOWDP(column) ! H2OSOI_ICE_update(layer,column) = H2OSOI_ICE(layer,column)+GainH2O_ICE(layer,column) ! H2OSOI_LIQ_update(layer,column) = H2OSOI_LIQ(layer,column)+GainH2O_LIQ(layer,column) ! ivar state ! 1 frac_sno ! 2 H2OSNO ! 3 H2OSOI_LIQ ! 4 H2OSOI_ICE ! 5 SNOWDP ! 6 DZSNO !!========================================================================= !--------------------------------------------------------------------------------------------------------kyh05072014 !--------------------------------------------------------------------------------------------------------kyh11142014 allocate(snlsno(Ncolumn)) call nc_check(nf90_inq_varid(ncFileID,'SNLSNO', VarID), 'sv_to_restart_file', 'inq_varid SNLSNO') call nc_check(nf90_get_var(ncFileID, VarID, snlsno), 'sv_to_restart_file', 'get_var SNLSNO') !--------------------------------------------------------------------------------------------------------kyh11142014 UPDATE : do ivar=1, nfields varname = trim(progvar(ivar)%varname) string2 = trim(filename)//' '//trim(varname) if ( .not. progvar(ivar)%update ) then write(string1,*)'intentionally not updating '//trim(string2) write(string3,*)'as per namelist control in model_nml:clm_variables' call error_handler(E_MSG, 'sv_to_restart_file', string1, text2=string3) cycle UPDATE endif if (trim(varname) == 'ZWT') then ! ZWT is calculated from WA and WT ... so we have to update those ! CLM variables based on the new ZWT from the assimilation. ! Simply updating ZWT will have no effect because upon restart ! CLM will calculate ZWT given the same old WA and WT. call update_water_table_depth( ivar, state_vector, ncFileID, filename, dart_time) cycle UPDATE endif ! Ensure netCDF variable is conformable with progvar quantity. ! The TIME and Copy dimensions are intentionally not queried ! by looping over the dimensions stored in the progvar type. call nc_check(nf90_inq_varid(ncFileID, varname, VarID), & 'sv_to_restart_file', 'inq_varid '//trim(string2)) call nc_check(nf90_inquire_variable(ncFileID,VarID,dimids=dimIDs,ndims=ncNdims), & 'sv_to_restart_file', 'inquire '//trim(string2)) DimCheck : do i = 1,progvar(ivar)%numdims write(string1,'(''inquire dimension'',i2,A)') i,trim(string2) call nc_check(nf90_inquire_dimension(ncFileID, dimIDs(i), len=dimlen), & 'sv_to_restart_file', string1) if ( dimlen /= progvar(ivar)%dimlens(i) ) then write(string1,*) trim(string2),' dim/dimlen ',i,dimlen,' not ',progvar(ivar)%dimlens(i) write(string2,*)' but it should be.' call error_handler(E_ERR, 'sv_to_restart_file', string1, & source, revision, revdate, text2=string2) endif enddo DimCheck ! When called with a 4th argument, vector_to_prog_var() replaces the DART ! missing code with the value in the corresponding variable in the netCDF file. ! Any clamping to physically meaningful values occurrs in vector_to_prog_var. if (progvar(ivar)%numdims == 1) then ni = progvar(ivar)%dimlens(1) !------------------------------------------------------------------------------original !allocate(data_1d_array(ni)) !call vector_to_prog_var(state_vector, ivar, data_1d_array, ncFileID) ! !call nc_check(nf90_put_var(ncFileID, VarID, data_1d_array), & ! 'sv_to_restart_file', 'put_var '//trim(varname)) !deallocate(data_1d_array) !------------------------------------------------------------------------------original !-------------------------------------------------------------------------------kyh11142014 allocate(data_1d_array(ni)) allocate(data_1d_array_prior(ni)) call vector_to_prog_var(state_vector, ivar, data_1d_array, data_1d_array_prior, ncFileID) if (progvar(ivar)%varname == 'H2OSNO') then write(*,*)'reading H2OSNO' allocate(h2osno_pr(ni)) allocate(h2osno_po(ni)) h2osno_po = data_1d_array h2osno_pr = data_1d_array_prior deallocate(data_1d_array) deallocate(data_1d_array_prior) elseif (progvar(ivar)%varname == 'SNOWDP') then write(*,*)'reading SNOWDP' allocate(snowdp_pr(ni)) allocate(snowdp_po(ni)) allocate(gain_snowdp(ni)) snowdp_po = data_1d_array snowdp_pr = data_1d_array_prior deallocate(data_1d_array) deallocate(data_1d_array_prior) else call nc_check(nf90_put_var(ncFileID, VarID, data_1d_array), & 'sv_to_restart_file', 'put_var '//trim(varname)) deallocate(data_1d_array) deallocate(data_1d_array_prior) endif !-------------------------------------------------------------------------------kyh11142014 elseif (progvar(ivar)%numdims == 2) then ni = progvar(ivar)%dimlens(1) nj = progvar(ivar)%dimlens(2) !-------------------------------------------------------------------------------original !allocate(data_2d_array(ni, nj)) !call vector_to_prog_var(state_vector, ivar, data_2d_array, ncFileID) ! !call nc_check(nf90_put_var(ncFileID, VarID, data_2d_array), & ! 'sv_to_restart_file', 'put_var '//trim(varname)) !deallocate(data_2d_array) !-------------------------------------------------------------------------------original !-------------------------------------------------------------------------------kyh11142014 allocate(data_2d_array(ni, nj)) allocate(data_2d_array_prior(ni, nj)) call vector_to_prog_var(state_vector, ivar, data_2d_array, data_2d_array_prior, ncFileID) if (progvar(ivar)%varname == "DZSNO") then write(*,*)'reading DZSNO' allocate(dzsno_pr(ni, nj)) allocate(dzsno_po(ni, nj)) !Temp variables allocated here allocate(gain_h2osno_l(ni,nj)) ! the SWE gained in each layer (total SWE gained *weight) allocate(gain_dzsno(ni,nj)) ! the gained snow thickness of each layer allocate(gain_h2oliq_l(ni,nj)) ! the gained snow water of each layer allocate(gain_h2oice_l(ni,nj)) ! the gained snow ice of each layer dzsno_pr = data_2d_array_prior deallocate(data_2d_array) deallocate(data_2d_array_prior) elseif (progvar(ivar)%varname == "H2OSOI_LIQ") then write(*,*)'reading H2OSOI_LIQ' allocate(h2oliq_pr(ni, nj)) allocate(h2oliq_po(ni, nj)) h2oliq_pr = data_2d_array_prior h2oliq_po = data_2d_array deallocate(data_2d_array) deallocate(data_2d_array_prior) elseif(progvar(ivar)%varname == "H2OSOI_ICE") then write(*,*)'reading H2OSOI_ICE' allocate(h2oice_pr(ni, nj)) allocate(h2oice_po(ni, nj)) h2oice_pr = data_2d_array_prior h2oice_po = data_2d_array deallocate(data_2d_array) deallocate(data_2d_array_prior) else call nc_check(nf90_put_var(ncFileID, VarID, data_2d_array), & 'sv_to_restart_file', 'put_var '//trim(varname)) deallocate(data_2d_array) deallocate(data_2d_array_prior) endif !-------------------------------------------------------------------------------kyh11142014 else write(string1, *) 'no support for data array of dimension ', ncNdims call error_handler(E_ERR,'sv_to_restart_file', string1, & source,revision,revdate) endif ! TJH FIXME ... this works perfectly if it were not for a bug in netCDF. ! When they fix the bug, this will be a useful thing to restore. ! Make note that the variable has been updated by DART ! call nc_check(nf90_Redef(ncFileID),'sv_to_restart_file', 'redef '//trim(filename)) ! call nc_check(nf90_put_att(ncFileID, VarID,'DART','variable modified by DART'),& ! 'sv_to_restart_file', 'modified '//trim(varname)) ! call nc_check(nf90_enddef(ncfileID),'sv_to_restart_file','state enddef '//trim(filename)) enddo UPDATE !-------------------------------------------------------------------------------------kyh11142014 !h2oliq_po = h2oliq_pr !h2oice_po = h2oice_pr dzsno_po = dzsno_pr snowdp_po = snowdp_pr wt_liq_l = 0 snowden = 0 !We have read all the temp vars needed. Now we need to adjust the variables to be consistant with !the updated H2OSNO !Now we have more than just H2OSNO to update gain_snowdp = 0 reduce_factor = 0.1_r8 do j = 1,nj if (h2osno_po(j)<0) then write(*,*)'negative value found: column ',j if(snlsno(j)==0) h2osno_po(j) = reduce_factor*h2osno_pr(j) do i=1,-snlsno(j) h2oliq_po(i,j)= reduce_factor*h2oliq_pr(i,j) h2oice_po(i,j)= reduce_factor*h2oice_pr(i,j) dzsno_po(i,j) = reduce_factor*dzsno_pr(i,j) !FIX ME. The value chosen is kind of random. end do snowdp_po(j) = sum(dzsno_po(:,j)) h2osno_po(j) = sum(h2oice_po(:,j))+sum(h2oliq_po(:,j)) else if(snlsno(j)<0)then ! write(*,*) 'snlsno(j)',snlsno(j) do c = 1, -snlsno(j) i = 5 - c +1 ! write(*,*)'================column',j,'==============================' ! write(*,*)'h2oliq_pr(i,j) in layer',i,':',h2oliq_pr(i,j) ! write(*,*)'h2oice_pr(i,j) in layer',i,':',h2oice_pr(i,j) ! write(*,*)'dzsno_pr(i,j) in layer',i,':',dzsno_pr(i,j) ! write(*,*)'h2osno(j)',h2osno_pr(j) wt_liq_l = h2oliq_pr(i,j)/(h2oliq_pr(i,j) + h2oice_pr(i,j)) ! write(*,*)'wt_liq_l',wt_liq_l snowden = (h2oliq_pr(i,j) + h2oice_pr(i,j))/dzsno_pr(i,j) ! write(*,*)'snoden',snowden if(snowden == NaN) snowden = 0.0_r8 wt_swe_l = (h2oliq_pr(i,j) + h2oice_pr(i,j))/h2osno_pr(j) if(wt_swe_l == NaN) wt_swe_l = 0.0_r8 ! write(*,*)'wt_swe_l',wt_swe_l gain_h2osno_l(i,j) = (h2osno_po(j) - h2osno_pr(j))*wt_swe_l gain_h2oliq_l(i,j) = gain_h2osno_l(i,j)*wt_liq_l gain_h2oice_l(i,j) = gain_h2osno_l(i,j)*(1 - wt_liq_l) gain_dzsno(i,j) = gain_h2osno_l(i,j)/snowden if(gain_dzsno(i,j) == NaN) gain_dzsno(i,j) =0.0_r8 h2oliq_po(i,j) = h2oliq_pr(i,j) + gain_h2oliq_l(i,j) h2oice_po(i,j) = h2oice_pr(i,j) + gain_h2oice_l(i,j) dzsno_po(i,j) = dzsno_pr(i,j) + gain_dzsno(i,j) end do gain_snowdp(j) = sum(gain_dzsno(:,j)) snowdp_po(j) = snowdp_pr(j) + gain_snowdp(j) endif endif enddo j = 6139 write(*,*)'snowdp_po:',snowdp_po(j) write(*,*)'gain_snowdp:',gain_snowdp(j) write(*,*)'snowdp_pr:',snowdp_pr(j) write(*,*)'snlsno(j):',snlsno(j) write(*,*)'h2osno_pr',h2osno_pr(j) write(*,*)'h2osno_po',h2osno_po(j) if(snlsno(j)<0)then ! write(*,*) 'snlsno(j)',snlsno(j) do c = 1, -snlsno(j) i = 5 - c +1 write(*,*)'h2oliq_po:',h2oliq_po(i,j) write(*,*)'h2oliq_pr:',h2oliq_pr(i,j) end do endif deallocate(gain_h2osno_l) deallocate(gain_h2oliq_l) deallocate(gain_h2oice_l) deallocate(gain_dzsno) deallocate(gain_snowdp) deallocate(snowdp_pr) deallocate(h2osno_pr) deallocate(h2oliq_pr) deallocate(h2oice_pr) deallocate(dzsno_pr) !UPDATE call nc_check(nf90_inq_varid(ncFileID,'H2OSNO', VarID), 'sv_to_restart_file', 'inq_varid H2OSNO') call nc_check(nf90_put_var(ncFileID, VarID, h2osno_po), & 'sv_to_restart_file', 'put_var '//'H2OSNO') deallocate(h2osno_po) call nc_check(nf90_inq_varid(ncFileID,'SNOWDP', VarID), 'sv_to_restart_file', 'inq_varid SNOWDP') call nc_check(nf90_put_var(ncFileID, VarID, snowdp_po), & 'sv_to_restart_file', 'put_var '//'SNOWDP') deallocate(snowdp_po) call nc_check(nf90_inq_varid(ncFileID,'DZSNO', VarID), 'sv_to_restart_file', 'inq_varid DZSNO') call nc_check(nf90_put_var(ncFileID, VarID, dzsno_po), & 'sv_to_restart_file', 'put_var '//'DZSNO') deallocate(dzsno_po) call nc_check(nf90_inq_varid(ncFileID,'H2OSOI_LIQ', VarID), 'sv_to_restart_file', 'inq_varid H2OSOI_LIQ') call nc_check(nf90_put_var(ncFileID, VarID, h2oliq_po), & 'sv_to_restart_file', 'put_var '//'H2OSOI_LIQ') deallocate(h2oliq_po) call nc_check(nf90_inq_varid(ncFileID,'H2OSOI_ICE', VarID), 'sv_to_restart_file', 'inq_varid H2OSOI_ICE') call nc_check(nf90_put_var(ncFileID, VarID, h2oice_po), & 'sv_to_restart_file', 'put_var '//'H2OSOI_ICE') deallocate(h2oice_po) !-------------------------------------------------------------------------------------kyh11142014 call nc_check(nf90_close(ncFileID),'sv_to_restart_file','close '//trim(filename)) ncFileID = 0 end subroutine sv_to_restart_file function get_ncols_in_gridcell(ilon,ilat) integer, intent(in) :: ilon, ilat integer :: get_ncols_in_gridcell if ( .not. module_initialized ) call static_init_model get_ncols_in_gridcell = gridCellInfo(ilon,ilat)%ncols end function get_ncols_in_gridcell subroutine get_colids_in_gridcell(ilon,ilat,colids) integer, intent(in) :: ilon, ilat integer, dimension(:), intent(out) :: colids if ( .not. module_initialized ) call static_init_model if (size(colids) == size(gridCellInfo(ilon,ilat)%columnids) ) then colids(:) = gridCellInfo(ilon,ilat)%columnids else write(string1,*)'lengths are wrong',size(colids) call error_handler(E_ERR,'get_colids_in_gridcell',string1,source,revision,revdate) endif end subroutine get_colids_in_gridcell !================================================================== ! The remaining interfaces come last !================================================================== subroutine model_interpolate(x, location, obs_kind, interp_val, istatus, optionals) ! PURPOSE: ! ! For a given lat, lon, and height, interpolate the correct state value ! to that location for the filter from the clm state vectors ! ! Reconstructing the vertical profile of the gridcell is complicated. ! Each land unit/column can have a different number of vertical levels. ! Do I just try to recontruct the whole profile and then interpolate? ! Impossible to know which elements are 'above' and 'below' without ! finding all the elements in the first place. The vertical information ! is in the levels() array for each state vector component. ! Passed variables real(r8), intent(in) :: x(:) type(location_type), intent(in) :: location integer, intent(in) :: obs_kind real(r8), intent(out) :: interp_val integer, intent(out) :: istatus real(r8), dimension(:), optional, intent(in) :: optionals ! Local storage real(r8), dimension(LocationDims) :: loc_array real(r8) :: llon, llat, lheight real(r8) :: interp_val_2 integer :: istatus_2 if ( .not. module_initialized ) call static_init_model ! Let's assume failure. Set return val to missing, then the code can ! just set istatus to something indicating why it failed, and return. ! If the interpolation is good, the interp_val will be set to the ! good value, and the last line here sets istatus to 0. ! make any error codes set here be in the 10s interp_val = MISSING_R8 ! the DART bad value flag interp_val_2 = MISSING_R8 ! the DART bad value flag istatus = 99 ! unknown error ! Get the individual locations values loc_array = get_location(location) llon = loc_array(1) llat = loc_array(2) lheight = loc_array(3) if ((debug > 6) .and. do_output()) print *, 'requesting interpolation at ', llon, llat, lheight ! FIXME may be better to check the %maxlevels and kick the interpolation to the ! appropriate routine based on that ... or check the dimnames for the ! vertical coordinate ... if (obs_kind == KIND_SOIL_TEMPERATURE) then call get_grid_vertval(x, location, 'T_SOISNO', interp_val, istatus ) elseif (obs_kind == KIND_SOIL_MOISTURE) then ! TJH FIXME - actually ROLAND FIXME ! This is terrible ... the COSMOS operator wants m3/m3 ... CLM is kg/m2 call get_grid_vertval(x, location, 'H2OSOI_LIQ',interp_val , istatus ) call get_grid_vertval(x, location, 'H2OSOI_ICE',interp_val_2, istatus_2 ) if ((istatus == 0) .and. (istatus_2 == 0)) then interp_val = interp_val + interp_val_2 else interp_val = MISSING_R8 istatus = 6 endif elseif (obs_kind == KIND_BRIGHTNESS_TEMPERATURE ) then if (present(optionals)) then call get_brightness_temperature(model_time, location, optionals, interp_val, istatus) else write(string1, '(''Tb obs at lon,lat ('',f12.6,'','',f12.6,'') has no metadata.'')') & llon, llat write(string2,*)'cannot call model_interpolate() for Tb without metadata argument.' call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate,text2=string2) endif elseif (obs_kind == KIND_LIQUID_WATER ) then call get_grid_vertval(x, location, 'H2OSOI_LIQ',interp_val, istatus ) elseif (obs_kind == KIND_ICE ) then call get_grid_vertval(x, location, 'H2OSOI_ICE',interp_val, istatus ) elseif (obs_kind == KIND_SNOWCOVER_FRAC ) then call compute_gridcell_value(x, location, 'frac_sno', interp_val, istatus) elseif (obs_kind == KIND_LEAF_CARBON ) then call compute_gridcell_value(x, location, 'leafc', interp_val, istatus) elseif (obs_kind == KIND_WATER_TABLE_DEPTH ) then call compute_gridcell_value(x, location, 'ZWT', interp_val, istatus) elseif (obs_kind == KIND_SNOW_THICKNESS ) then write(string1,*)'model_interpolate for DZSNO not written yet.' call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate) istatus = 5 elseif ((obs_kind == KIND_GEOPOTENTIAL_HEIGHT) .and. vert_is_level(location)) then if (nint(lheight) > nlevgrnd) then interp_val = MISSING_R8 istatus = 1 else interp_val = LEVGRND(nint(lheight)) istatus = 0 endif else write(string1,*)'model_interpolate not written for (integer) kind ',obs_kind call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate) istatus = 5 endif if ((debug > 6) .and. do_output()) write(*,*)'interp_val ',interp_val end subroutine model_interpolate !------------------------------------------------------------------ subroutine compute_gridcell_value(x, location, varstring, interp_val, istatus) ! ! Each gridcell may contain values for several land units, each land unit may contain ! several columns, each column may contain several pft's. BUT this routine never ! aggregates across multiple pft's. So, each gridcell value ! is an area-weighted value of an unknown number of column-based quantities. ! Passed variables real(r8), intent(in) :: x(:) ! state vector type(location_type), intent(in) :: location ! location somewhere in a grid cell character(len=*), intent(in) :: varstring ! frac_sno, leafc real(r8), intent(out) :: interp_val ! area-weighted result integer, intent(out) :: istatus ! error code (0 == good) ! Local storage integer :: ivar, index1, indexN, indexi, counter integer :: gridloni,gridlatj real(r8) :: loc_lat, loc_lon real(r8) :: total, total_area real(r8), dimension(1) :: loninds,latinds real(r8), dimension(LocationDims) :: loc if ( .not. module_initialized ) call static_init_model ! Let's assume failure. Set return val to missing, then the code can ! just set istatus to something indicating why it failed, and return. ! If the interpolation is good, the interp_val will be set to the ! good value, and the last line here sets istatus to 0. ! make any error codes set here be in the 10s interp_val = MISSING_R8 ! the DART bad value flag istatus = 99 ! unknown error loc = get_location(location) ! loc is in DEGREES loc_lon = loc(1) loc_lat = loc(2) ! determine the portion of interest of the state vector ivar = findVarIndex(varstring, 'compute_gridcell_value') index1 = progvar(ivar)%index1 ! in the DART state vector, start looking here indexN = progvar(ivar)%indexN ! in the DART state vector, stop looking here ! BOMBPROOFING - check for a vertical dimension for this variable if (progvar(ivar)%maxlevels > 1) then 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, & source, revision, revdate, text2=string2) endif ! determine the grid cell for the location latinds = minloc(abs(LAT - loc_lat)) ! these return 'arrays' ... loninds = minloc(abs(LON - loc_lon)) ! these return 'arrays' ... gridlatj = latinds(1) gridloni = loninds(1) if ((debug > 5) .and. do_output()) then write(*,*)'compute_gridcell_value:targetlon, lon, lon index is ',& loc_lon,LON(gridloni),gridloni write(*,*)'compute_gridcell_value:targetlat, lat, lat index is ',& loc_lat,LAT(gridlatj),gridlatj endif ! If there is no vertical component, the problem is greatly simplified. ! Simply area-weight an average of all pieces in the grid cell. ! FIXME ... this is the loop that can exploit the knowledge of what ! columnids or pftids are needed for any particular gridcell. ! gridCellInfo%pftids, gridCellInfo%columnids counter = 0 total = 0.0_r8 ! temp storage for state vector total_area = 0.0_r8 ! temp storage for area ELEMENTS : do indexi = index1, indexN if ( lonixy(indexi) /= gridloni ) cycle ELEMENTS if ( latjxy(indexi) /= gridlatj ) cycle ELEMENTS if ( x(indexi) == MISSING_R8) cycle ELEMENTS if ( landarea(indexi) == 0.0_r8 ) cycle ELEMENTS counter = counter + 1 total = total + x(indexi)*landarea(indexi) total_area = total_area + landarea(indexi) if ((debug > 5) .and. do_output()) then write(*,*) write(*,*)'gridcell location match',counter,'at statevector index',indexi write(*,*)'statevector value is (',x(indexi),')' write(*,*)'area is (',landarea(indexi),')' write(*,*)'LON index is (',lonixy(indexi),')' write(*,*)'LAT index is (',latjxy(indexi),')' write(*,*)'closest LON is (',LON(gridloni),')' write(*,*)'closest LAT is (',LAT(gridlatj),')' write(*,*)'closest lev is (',levels(indexi),')' endif enddo ELEMENTS if (total_area /= 0.0_r8) then ! All good. interp_val = total/total_area istatus = 0 else if ((debug > 4) .and. do_output()) then write(string1, *)'Variable '//trim(varstring)//' had no viable data' write(string2, *)'at gridcell ilon/jlat = (',gridloni,',',gridlatj,')' write(string3, *)'obs lon/lat = (',loc_lon,',',loc_lat,')' call error_handler(E_MSG,'compute_gridcell_value', string1, & text2=string2,text3=string3) endif endif ! Print more information for the really curious if ((debug > 5) .and. do_output()) then write(string1,*)'counter, total, total_area', counter, total, total_area write(string2,*)'interp_val, istatus', interp_val, istatus call error_handler(E_MSG,'compute_gridcell_value', string1, text2=string2) endif end subroutine compute_gridcell_value !------------------------------------------------------------------ subroutine get_grid_vertval(x, location, varstring, interp_val, istatus) ! ! Calculate the expected vertical value for the gridcell. ! Each gridcell value is an area-weighted value of an unknown number of ! column-based quantities. ! Passed variables real(r8), intent(in) :: x(:) ! state vector type(location_type), intent(in) :: location ! location somewhere in a grid cell character(len=*), intent(in) :: varstring ! T_SOISNO, H2OSOI_LIQ, H2OSOI_ICE real(r8), intent(out) :: interp_val ! area-weighted result integer, intent(out) :: istatus ! error code (0 == good) ! Local storage integer :: ivar, index1, indexN, indexi, counter1, counter2 integer :: gridloni,gridlatj real(r8), dimension(LocationDims) :: loc real(r8) :: loc_lat, loc_lon, loc_lev real(r8) :: value_below, value_above, total_area real(r8) :: depthbelow, depthabove real(r8) :: topwght, botwght real(r8), dimension(1) :: loninds,latinds real(r8), allocatable, dimension(:) :: above, below real(r8), allocatable, dimension(:,:) :: myarea if ( .not. module_initialized ) call static_init_model ! Let's assume failure. Set return val to missing, then the code can ! just set istatus to something indicating why it failed, and return. ! If the interpolation is good, the interp_val will be set to the ! good value, and the last line here sets istatus to 0. ! make any error codes set here be in the 10s interp_val = MISSING_R8 ! the DART bad value flag istatus = 99 ! unknown error loc = get_location(location) ! loc is in DEGREES loc_lon = loc(1) loc_lat = loc(2) loc_lev = loc(3) if ( loc_lev < 0.0_r8 ) then write(string1,*)'Cannot support above-ground vertical interpolation.' write(string2,*)'requested a value at a depth of ',loc_lev write(string3,*)'CLM has negative depths to indicate above-ground values.' call error_handler(E_ERR,'get_grid_vertval', string1, & source, revision, revdate, text2=string2, text3=string3) endif ! determine the portion of interest of the state vector ivar = findVarIndex(varstring, 'get_grid_vertval') index1 = progvar(ivar)%index1 ! in the DART state vector, start looking here indexN = progvar(ivar)%indexN ! in the DART state vector, stop looking here ! BOMBPROOFING - check for a vertical dimension for this variable if (progvar(ivar)%maxlevels < 2) then write(string1, *)'Variable '//trim(varstring)//' should not use this routine.' write(string2, *)'use compute_gridcell_value() instead.' call error_handler(E_ERR,'get_grid_vertval', string1, & source, revision, revdate, text2=string2) endif ! determine the grid cell for the location latinds = minloc(abs(LAT - loc_lat)) ! these return 'arrays' ... loninds = minloc(abs(LON - loc_lon)) ! these return 'arrays' ... gridlatj = latinds(1) gridloni = loninds(1) if ((debug > 4) .and. do_output()) then write(*,*)'get_grid_vertval:targetlon, lon, lon index, level is ', & loc_lon,LON(gridloni),gridloni,loc_lev write(*,*)'get_grid_vertval:targetlat, lat, lat index, level is ', & loc_lat,LAT(gridlatj),gridlatj,loc_lev endif ! Determine the level 'above' and 'below' the desired vertical ! The above-ground 'depths' are calculated from ZISNO and are negative. ! The 'depths' are all positive numbers, increasingly positive is deeper. ! The variables currently supported use the subsurface definitions in ! the module variable LEVNGRND. if (loc_lev <= LEVGRND(1)) then ! the top level is so close to the surface depthabove = LEVGRND(1) ! just use the top level depthbelow = LEVGRND(1) elseif (loc_lev >= maxval(LEVGRND)) then ! at depth, however ... do we depthabove = maxval(LEVGRND) ! fail or just use the deepest depthbelow = maxval(LEVGRND) ! I am using the deepest. else LAYERS : do indexi = 2,size(LEVGRND) if (loc_lev < LEVGRND(indexi)) then depthabove = LEVGRND(indexi-1) depthbelow = LEVGRND(indexi ) exit LAYERS endif enddo LAYERS endif if ((debug > 4) .and. do_output()) then write(*,*)'get_grid_vertval:depthbelow ',depthbelow,'>= loc_lev', & loc_lev,'>= depthabove',depthabove endif ! Determine how many elements can contribute to the gridcell value. ! There are multiple column-based contributors, each column has a ! separate area-based weight. There are multiple levels. ! I believe I have to keep track of all of them to sort out how to ! calculate the gridcell value at a particular depth. counter1 = 0 counter2 = 0 GRIDCELL : do indexi = index1, indexN if ( lonixy(indexi) /= gridloni ) cycle GRIDCELL if ( latjxy(indexi) /= gridlatj ) cycle GRIDCELL if ( x(indexi) == MISSING_R8) cycle GRIDCELL if (levels(indexi) == depthabove) counter1 = counter1 + 1 if (levels(indexi) == depthbelow) counter2 = counter2 + 1 enddo GRIDCELL if ( (counter1+counter2) == 0 ) then if ((debug > 0) .and. do_output()) then write(string1, *)'statevector variable '//trim(varstring)//' had no viable data' write(string2, *)'at gridcell lon/lat = (',gridloni,',',gridlatj,')' write(string3, *)'obs lon/lat/lev (',loc_lon,',',loc_lat,',',loc_lev,')' call error_handler(E_MSG,'get_grid_vertval', string1, & text2=string2,text3=string3) endif return endif allocate(above(counter1),below(counter2),myarea(max(counter1,counter2),2)) above = MISSING_R8 below = MISSING_R8 myarea = 0.0_r8 counter1 = 0 counter2 = 0 ELEMENTS : do indexi = index1, indexN if ( lonixy(indexi) /= gridloni ) cycle ELEMENTS if ( latjxy(indexi) /= gridlatj ) cycle ELEMENTS if ( x(indexi) == MISSING_R8) cycle ELEMENTS ! write(*,*)'level ',indexi,' is ',levels(indexi),' location depth is ',loc_lev if (levels(indexi) == depthabove) then counter1 = counter1 + 1 above( counter1) = x(indexi) myarea(counter1,1) = landarea(indexi) endif if (levels(indexi) == depthbelow) then counter2 = counter2 + 1 below( counter2) = x(indexi) myarea(counter2,2) = landarea(indexi) endif if ((levels(indexi) /= depthabove) .and. & (levels(indexi) /= depthbelow)) then cycle ELEMENTS endif if ((debug > 4) .and. do_output()) then write(*,*) write(*,*)'gridcell location match at statevector index',indexi write(*,*)'statevector value is (',x(indexi),')' write(*,*)'area is (',landarea(indexi),')' write(*,*)'LON index is (',lonixy(indexi),')' write(*,*)'LAT index is (',latjxy(indexi),')' write(*,*)'gridcell LON is (',LON(gridloni),')' write(*,*)'gridcell LAT is (',LAT(gridlatj),')' write(*,*)'depth is (',levels(indexi),')' endif enddo ELEMENTS ! could arise if the above or below was 'missing' ... but the mate was not. if ( counter1 /= counter2 ) then write(string1, *)'Variable '//trim(varstring)//' has peculiar interpolation problems.' write(string2, *)'uneven number of values "above" and "below"' write(string3, *)'counter1 == ',counter1,' /= ',counter2,' == counter2' call error_handler(E_MSG,'get_grid_vertval', string1, & text2=string2,text3=string3) return endif ! Determine the value for the level above the depth of interest. total_area = sum(myarea(1:counter1,1)) if ( total_area /= 0.0_r8 ) then ! normalize the area-based weights myarea(1:counter1,1) = myarea(1:counter1,1) / total_area value_above = sum(above(1:counter1) * myarea(1:counter1,1)) else write(string1, *)'Variable '//trim(varstring)//' had no viable data above' write(string2, *)'at gridcell lon/lat/lev = (',gridloni,',',gridlatj,',',depthabove,')' write(string3, *)'obs lon/lat/lev (',loc_lon,',',loc_lat,',',loc_lev,')' call error_handler(E_ERR,'get_grid_vertval', string1, & source, revision, revdate, text2=string2,text3=string3) endif ! Determine the value for the level below the depth of interest. total_area = sum(myarea(1:counter2,2)) if ( total_area /= 0.0_r8 ) then ! normalize the area-based weights myarea(1:counter2,2) = myarea(1:counter2,2) / total_area value_below = sum(below(1:counter2) * myarea(1:counter2,2)) else write(string1, *)'Variable '//trim(varstring)//' had no viable data below' write(string2, *)'at gridcell lon/lat/lev = (',gridloni,',',gridlatj,',',depthbelow,')' write(string3, *)'obs lon/lat/lev (',loc_lon,',',loc_lat,',',loc_lev,')' call error_handler(E_ERR,'get_grid_vertval', string1, & source, revision, revdate, text2=string2,text3=string3) endif if (depthbelow == depthabove) then topwght = 1.0_r8 botwght = 0.0_r8 else topwght = (depthbelow - loc_lev) / (depthbelow - depthabove) botwght = (loc_lev - depthabove) / (depthbelow - depthabove) endif interp_val = value_above*topwght + value_below*botwght istatus = 0 deallocate(above, below, myarea) end subroutine get_grid_vertval !------------------------------------------------------------------ !subroutine vector_to_1d_prog_var(x, ivar, data_1d_array, ncid) !original subroutine vector_to_1d_prog_var(x, ivar, data_1d_array, data_1d_array_prior, ncid) !kyh11142014 !------------------------------------------------------------------ ! convert the values from a 1d array, starting at an offset, into a 1d array. ! ! If the optional argument (ncid) is specified, some additional ! processing takes place. The variable in the netcdf is read. ! This must be the same shape as the intended output array. ! Anywhere the DART MISSING code is encountered in the input array, ! the corresponding (i.e. original) value from the netCDF file is ! used. real(r8), dimension(:), intent(in) :: x integer, intent(in) :: ivar real(r8), dimension(:), intent(out) :: data_1d_array real(r8), dimension(:), intent(out) :: data_1d_array_prior !kyh11142014 integer, OPTIONAL, intent(in) :: ncid integer :: i,ii, VarID real(r8), allocatable, dimension(:) :: org_array ! unpack the right part of the DART state vector into a 1D array. ii = progvar(ivar)%index1 do i = 1, progvar(ivar)%dimlens(1) data_1d_array(i) = x(ii) ii = ii + 1 enddo ii = ii - 1 if ( ii /= progvar(ivar)%indexN ) then write(string1, *)'Variable '//trim(progvar(ivar)%varname)//' filled wrong.' write(string2, *)'Should have ended at ',progvar(ivar)%indexN,' actually ended at ',ii call error_handler(E_ERR,'vector_to_1d_prog_var', string1, & source, revision, revdate, text2=string2) endif ! Apply the min/max values, if applicable ! This should only be true when converting to a variable that will ! be reinserted into the CLM restart file. This is indicated ! by the presence of the ncid variable. if (present(ncid)) then if ((progvar(ivar)%rangeRestricted == BOUNDED_ABOVE ) .or. & (progvar(ivar)%rangeRestricted == BOUNDED_BOTH )) then where ((data_1d_array /= MISSING_R8) .and. & (data_1d_array > progvar(ivar)%maxvalue)) & data_1d_array = progvar(ivar)%maxvalue endif if ((progvar(ivar)%rangeRestricted == BOUNDED_BELOW ) .or. & (progvar(ivar)%rangeRestricted == BOUNDED_BOTH )) then where ((data_1d_array /= MISSING_R8) .and. & (data_1d_array < progvar(ivar)%minvalue)) & data_1d_array = progvar(ivar)%minvalue endif ! Replace the DART fill value with the original value and apply any clamping. ! Get the 'original' variable from the netcdf file. allocate(org_array(size(data_1d_array))) call nc_check(nf90_inq_varid(ncid, progvar(ivar)%varname, VarID), & 'vector_to_1d_prog_var', 'inq_varid '//trim(progvar(ivar)%varname)) call nc_check(nf90_get_var(ncid, VarID, org_array), & 'vector_to_1d_prog_var', 'get_var '//trim(progvar(ivar)%varname)) ! restoring the indeterminate original values where(data_1d_array == MISSING_R8) data_1d_array = org_array data_1d_array_prior = org_array !kyh11142014 ! clamping the assimilated values to physically meaningful ranges. !---------------------------------------------------------------------------original !if (trim(progvar(ivar)%varname) == 'SNOWDP') & ! where((data_1d_array < 0.0_r8)) data_1d_array = org_array ! !if (trim(progvar(ivar)%varname) == 'H2OSNO') & ! where((data_1d_array <= 0.0_r8)) data_1d_array = org_array !---------------------------------------------------------------------------original !----------------------------------------------------------------------------------------------------------------------------------kyh01052014 (RTM parameters) if (trim(progvar(ivar)%varname) == 'SNOWDP') then where((data_1d_array < 0.0_r8)) data_1d_array = org_array elseif (trim(progvar(ivar)%varname) == 'H2OSNO') then where((data_1d_array <= 0.0_r8)) data_1d_array = org_array !---------------------------------------------------------------------kyh05062014 where((data_1d_array > 1000._r8)) data_1d_array = org_array !---------------------------------------------------------------------kyh05062014 !-------------------------------------------------------------------------------------------------kyh12182014 (RTM parameters) elseif (trim(progvar(ivar)%varname) == 'STICKINESS_C') then where((data_1d_array < 0.1_r8)) data_1d_array = org_array where((data_1d_array > 0.5_r8)) data_1d_array = org_array elseif (trim(progvar(ivar)%varname) == 'B_PRIME_C') then where((data_1d_array < 0.496_r8)) data_1d_array = org_array where((data_1d_array > 0.744_r8)) data_1d_array = org_array elseif (trim(progvar(ivar)%varname) == 'X_LAMBDA_C') then where((data_1d_array < -1.656_r8)) data_1d_array = org_array where((data_1d_array > -0.804_r8)) data_1d_array = org_array !endif !-------------------------------------------------------------------------------------------------kyh12182014 (RTM parameters) !-------------------------------------------------------------------------------------------------kyh03052016 elseif (trim(progvar(ivar)%varname) == 'SCATT_ALBEDO_C') then where((data_1d_array < 0.05_r8)) data_1d_array = org_array where((data_1d_array > 0.1_r8)) data_1d_array = org_array endif !-------------------------------------------------------------------------------------------------kyh03052016 !----------------------------------------------------------------------------------------------------------------------------------kyh01052014 (RTM parameters) deallocate(org_array) else if (progvar(ivar)%xtype == NF90_INT) then where(data_1d_array == MISSING_I) data_1d_array = progvar(ivar)%spvalINT elseif (progvar(ivar)%xtype == NF90_FLOAT) then where(data_1d_array == MISSING_R4) data_1d_array = progvar(ivar)%spvalR4 elseif (progvar(ivar)%xtype == NF90_DOUBLE) then where(data_1d_array == MISSING_R8) data_1d_array = progvar(ivar)%spvalR8 endif endif end subroutine vector_to_1d_prog_var !------------------------------------------------------------------ !subroutine vector_to_2d_prog_var(x, ivar, data_2d_array, ncid) !original subroutine vector_to_2d_prog_var(x, ivar, data_2d_array, data_2d_array_prior, ncid) !kyh11142014 !------------------------------------------------------------------ ! convert the values from a 1d array, starting at an offset, ! into a 2d array. ! real(r8), dimension(:), intent(in) :: x integer, intent(in) :: ivar real(r8), dimension(:,:), intent(out) :: data_2d_array real(r8), dimension(:,:), intent(out) :: data_2d_array_prior !kyh11142014 integer, OPTIONAL, intent(in) :: ncid integer :: i,j,ii, VarID real(r8), allocatable, dimension(:,:) :: org_array ! unpack the right part of the DART state vector into a 1D array. ii = progvar(ivar)%index1 do j = 1,progvar(ivar)%dimlens(2) do i = 1,progvar(ivar)%dimlens(1) data_2d_array(i,j) = x(ii) ii = ii + 1 enddo enddo ii = ii - 1 if ( ii /= progvar(ivar)%indexN ) then write(string1, *)'Variable '//trim(progvar(ivar)%varname)//' filled wrong.' write(string2, *)'Should have ended at ',progvar(ivar)%indexN,' actually ended at ',ii call error_handler(E_ERR,'vector_to_2d_prog_var', string1, & source, revision, revdate, text2=string2) endif ! Apply the min/max values, if applicable ! This should only be true when converting to a variable that will ! be reinserted into the CLM restart file. This is indicated ! by the presence of the ncid variable. if (present(ncid)) then if ((progvar(ivar)%rangeRestricted == BOUNDED_ABOVE ) .or. & (progvar(ivar)%rangeRestricted == BOUNDED_BOTH )) then where ((data_2d_array /= MISSING_R8) .and. & (data_2d_array > progvar(ivar)%maxvalue)) & data_2d_array = progvar(ivar)%maxvalue endif if ((progvar(ivar)%rangeRestricted == BOUNDED_BELOW ) .or. & (progvar(ivar)%rangeRestricted == BOUNDED_BOTH )) then where ((data_2d_array /= MISSING_R8) .and. & (data_2d_array < progvar(ivar)%minvalue)) & data_2d_array = progvar(ivar)%minvalue endif ! Replace the DART fill value with the original value and apply any clamping. ! Get the 'original' variable from the netcdf file if need be. allocate(org_array(size(data_2d_array,1),size(data_2d_array,2))) call nc_check(nf90_inq_varid(ncid, progvar(ivar)%varname, VarID), & 'vector_to_2d_prog_var', 'inq_varid '//trim(progvar(ivar)%varname)) call nc_check(nf90_get_var(ncid, VarID, org_array), & 'vector_to_2d_prog_var', 'get_var '//trim(progvar(ivar)%varname)) ! restoring the indeterminate original values where(data_2d_array == MISSING_R8 ) data_2d_array = org_array data_2d_array_prior = org_array !kyh11142014 ! clamping the assimilated values to physically meaningful ranges. if (trim(progvar(ivar)%varname) == 'DZSNO') then where((data_2d_array < 0.0_r8)) data_2d_array = org_array elseif (trim(progvar(ivar)%varname) == 'ZSNO') then where((data_2d_array < 0.0_r8)) data_2d_array = org_array 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 where((data_2d_array < 0.0_r8)) data_2d_array = org_array !------------------------------------------------------------------------kyh05052014 ! Currently, only update snow layers and soil surface layer, ! i.e., the 1th - 6th layer in levtot. ! So replace other layers' value with original value. !data_2d_array(8:20,:) = org_array(8:20,:) !data_2d_array(6:20,:) = org_array(6:20,:) !where((data_2d_array > 1.0_r8)) data_2d_array = org_array !------------------------------------------------------------------------kyh05052014 elseif (trim(progvar(ivar)%varname) == 'H2OSOI_ICE') then where((data_2d_array < 0.0_r8)) data_2d_array = org_array !------------------------------------------------------------------------kyh05052014 ! Currently, only update snow layers and soil surface layer, ! i.e., the 1th - 6th layer in levtot. ! So replace other layers' value with original value. !data_2d_array(8:20,:) = org_array(8:20,:) !data_2d_array(6:20,:) = org_array(6:20,:) !where((data_2d_array > 1.0_r8)) data_2d_array = org_array !------------------------------------------------------------------------kyh05052014 elseif (trim(progvar(ivar)%varname) == 'T_SOISNO') then where((data_2d_array < 0.0_r8)) data_2d_array = org_array !------------------------------------------------------------------------kyh05052014 ! Currently, only update snow layers and soil surface layer, ! i.e., the 1th - 6th layer in levtot. ! So replace other layers' value with original value. !data_2d_array(8:20,:) = org_array(8:20,:) !data_2d_array(6:20,:) = org_array(6:20,:) !------------------------------------------------------------------------kyh05052014 elseif (trim(progvar(ivar)%varname) == 'T_LAKE') then where((data_2d_array < 0.0_r8)) data_2d_array = org_array elseif (trim(progvar(ivar)%varname) == 'leafc') then where((data_2d_array < 0.0_r8)) data_2d_array = 0.0_r8 !---------------------------------------------------------------------------kyh05012014 elseif (trim(progvar(ivar)%varname) == 'snw_rds') then where((data_2d_array < 54.526_r8)) data_2d_array = org_array !microns where((data_2d_array > 1500._r8)) data_2d_array = org_array !microns !data_2d_array(1:4,:) = org_array(1:4,:) !kyh03122015 !---------------------------------------------------------------------------kyh05012014 endif deallocate(org_array) else if (progvar(ivar)%xtype == NF90_INT) then where(data_2d_array == MISSING_I) data_2d_array = progvar(ivar)%spvalINT elseif (progvar(ivar)%xtype == NF90_FLOAT) then where(data_2d_array == MISSING_R4) data_2d_array = progvar(ivar)%spvalR4 elseif (progvar(ivar)%xtype == NF90_DOUBLE) then where(data_2d_array == MISSING_R8) data_2d_array = progvar(ivar)%spvalR8 endif endif end subroutine vector_to_2d_prog_var !------------------------------------------------------------------ subroutine get_history_dims(ncid,fname,cstat,lon,lat,levgrnd,lonatm,latatm,lonrof,latrof) !------------------------------------------------------------------ ! ! Read the dimensions from the history netcdf file. ! ! The file name comes from module storage ... namelist. integer, intent(inout) :: ncid character(len=*), intent(in) :: fname character(len=*), intent(in) :: cstat ! how do you want to leave the netcdf file integer, intent(out) :: lon integer, intent(out) :: lat integer, intent(out) :: levgrnd integer, OPTIONAL, intent(out) :: lonatm, latatm integer, OPTIONAL, intent(out) :: lonrof, latrof integer :: dimid ! get the ball rolling ... if (ncid == 0) then ! we need to open it call nc_check(nf90_open(trim(fname), nf90_nowrite, ncid), & 'get_history_dims','open '//trim(fname)) endif ! The new SingleColumn (and unstructured grid) configurations ! do not have a 'lon' and 'lat' dimension. There is only 'lndgrid' if ( nf90_inq_dimid(ncid, 'lndgrid', dimid) == NF90_NOERR ) unstructured = .true. if (unstructured) then ! use the lndgrid dimension for both lon and lat call nc_check(nf90_inq_dimid(ncid, 'lndgrid', dimid), & 'get_history_dims','inq_dimid lndgrid '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=lon), & 'get_history_dims','inquire_dimension lndgrid '//trim(fname)) lat = lon else call nc_check(nf90_inq_dimid(ncid, 'lon', dimid), & 'get_history_dims','inq_dimid lon '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=lon), & 'get_history_dims','inquire_dimension lon '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'lat', dimid), & 'get_history_dims','inq_dimid lat '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=lat), & 'get_history_dims','inquire_dimension lat '//trim(fname)) endif call nc_check(nf90_inq_dimid(ncid, 'levgrnd', dimid), & 'get_history_dims','inq_dimid levgrnd '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=levgrnd), & 'get_history_dims','inquire_dimension levgrnd '//trim(fname)) if (present(lonatm)) then call nc_check(nf90_inq_dimid(ncid, 'lonatm', dimid), & 'get_history_dims','inq_dimid lonatm '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=lonatm), & 'get_history_dims','inquire_dimension lonatm '//trim(fname)) endif if (present(latatm)) then call nc_check(nf90_inq_dimid(ncid, 'latatm', dimid), & 'get_history_dims','inq_dimid latatm '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=latatm), & 'get_history_dims','inquire_dimension latatm '//trim(fname)) endif if (present(lonrof)) then call nc_check(nf90_inq_dimid(ncid, 'lonrof', dimid), & 'get_history_dims','inq_dimid lonrof '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=lonrof), & 'get_history_dims','inquire_dimension lonrof '//trim(fname)) endif if (present(latrof)) then call nc_check(nf90_inq_dimid(ncid, 'latrof', dimid), & 'get_history_dims','inq_dimid latrof '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=latrof), & 'get_history_dims','inquire_dimension latrof '//trim(fname)) endif if ((debug > 8) .and. do_output()) then write(logfileunit,*) write(logfileunit,*)'get_history_dims output follows:' write(logfileunit,*)'nlon = ',nlon write(logfileunit,*)'nlat = ',nlat write( * ,*) write( * ,*)'get_history_dims output follows:' write( * ,*)'nlon = ',nlon write( * ,*)'nlat = ',nlat if (present(lonatm)) write(logfileunit,*)'lonatm = ',lonatm if (present(latatm)) write(logfileunit,*)'latatm = ',latatm if (present(lonrof)) write(logfileunit,*)'lonrof = ',lonrof if (present(latrof)) write(logfileunit,*)'latrof = ',latrof if (present(lonatm)) write( * ,*)'lonatm = ',lonatm if (present(latatm)) write( * ,*)'latatm = ',latatm if (present(lonrof)) write( * ,*)'lonrof = ',lonrof if (present(latrof)) write( * ,*)'latrof = ',latrof endif if (cstat == 'close') then call nc_check(nf90_close(ncid),'get_history_dims','close '//trim(fname) ) ncid = 0 endif end subroutine get_history_dims !------------------------------------------------------------------ subroutine get_full_grid(ncid, fname, cstat) !------------------------------------------------------------------ ! ! Read the grid dimensions from the CLM history netcdf file. ! LON,LAT,AREA,LANDFRAC,LEVGRN all have module scope integer, intent(inout) :: ncid character(len=*), intent(in) :: fname character(len=*), intent(in) :: cstat ! Make sure the variables are the right size ... ! at some point in the future ... if (ncid == 0) then ! we need to open it call nc_check(nf90_open(trim(fname), nf90_nowrite, ncid), & 'get_full_grid','open '//trim(fname)) endif ! The lat/lon matrices in the history file have been masked by ! the land values such that the wet cells are 'missing' values. ! This makes it less than useful for us. Thankfully, the 1D ! lat/lon arrays have no such mask applied. We use these. call DART_get_var(ncid,'lon' ,LON) call DART_get_var(ncid,'lat' ,LAT) call DART_get_var(ncid,'levgrnd' ,LEVGRND) if (unstructured) then call DART_get_var(ncid,'area' ,AREA1D) call DART_get_var(ncid,'landfrac',LANDFRAC1D) where(AREA1D == MISSING_R8) AREA1D = 0.0_r8 where(LANDFRAC1D == MISSING_R8) LANDFRAC1D = 0.0_r8 else call DART_get_var(ncid,'area' ,AREA2D) call DART_get_var(ncid,'landfrac',LANDFRAC2D) where(AREA2D == MISSING_R8) AREA2D = 0.0_r8 where(LANDFRAC2D == MISSING_R8) LANDFRAC2D = 0.0_r8 endif ! just to make sure we are [0,360] and [-90,90] where (LON < 0.0_r8) LON = LON + 360.0_r8 where (LON > 360.0_r8) LON = LON - 360.0_r8 if (any(LON < 0.0_r8)) then write(string1,*)'longitudes in history file variable "lon" still negative.' call error_handler(E_ERR,'get_full_grid',string1,source,revision,revdate) endif where (LAT < -90.0_r8) LAT = -90.0_r8 where (LAT > 90.0_r8) LAT = 90.0_r8 if (cstat == 'close') then call nc_check(nf90_close(ncid),'get_full_grid','close '//trim(fname) ) ncid = 0 endif ! A little sanity check if ((debug > 7) .and. do_output()) then write(logfileunit,*) write(logfileunit,*)'history_file grid information as interpreted ...' write(logfileunit,*)'lon range ',minval(LON ),maxval(LON ) write(logfileunit,*)'lat range ',minval(LAT ),maxval(LAT ) write(logfileunit,*)'levgrnd range ',minval(LEVGRND ),maxval(LEVGRND ) write(logfileunit,*)'levgrnd is ',LEVGRND write( * ,*) write( * ,*)'history_file grid information as interpreted ...' write( * ,*)'lon range ',minval(LON ),maxval(LON ) write( * ,*)'lat range ',minval(LAT ),maxval(LAT ) write( * ,*)'levgrnd range ',minval(LEVGRND ),maxval(LEVGRND ) write( * ,*)'levgrnd is ',LEVGRND if (unstructured) then write(logfileunit,*)'area range ',minval(AREA1D ),maxval(AREA1D ) write(logfileunit,*)'landfrac range ',minval(LANDFRAC1D),maxval(LANDFRAC1D) write( * ,*)'area range ',minval(AREA1D ),maxval(AREA1D ) write( * ,*)'landfrac range ',minval(LANDFRAC1D),maxval(LANDFRAC1D) else write(logfileunit,*)'area range ',minval(AREA2D ),maxval(AREA2D ) write(logfileunit,*)'landfrac range ',minval(LANDFRAC2D),maxval(LANDFRAC2D) write( * ,*)'area range ',minval(AREA2D ),maxval(AREA2D ) write( * ,*)'landfrac range ',minval(LANDFRAC2D),maxval(LANDFRAC2D) endif endif return end subroutine get_full_grid !------------------------------------------------------------------ subroutine get_sparse_dims(ncid, fname, cstat) !------------------------------------------------------------------ ! ! Read the dimensions from the CLM restart netcdf file. integer, intent(inout) :: ncid character(len=*), intent(in) :: fname character(len=*), intent(in) :: cstat integer :: dimid, istatus, mylevgrnd if (ncid == 0) then call nc_check(nf90_open(trim(fname), nf90_nowrite, ncid), & 'get_sparse_dims','open '//trim(fname)) endif ! get dimid for 'gridcell' and then get value ... call nc_check(nf90_inq_dimid(ncid, 'gridcell', dimid), & 'get_sparse_dims','inq_dimid gridcell '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Ngridcell), & 'get_sparse_dims','inquire_dimension gridcell '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'landunit', dimid), & 'get_sparse_dims','inq_dimid landunit '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlandunit), & 'get_sparse_dims','inquire_dimension landunit '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'column', dimid), & 'get_sparse_dims','inq_dimid column '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Ncolumn), & 'get_sparse_dims','inquire_dimension column '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'pft', dimid), & 'get_sparse_dims','inq_dimid pft '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Npft), & 'get_sparse_dims','inquire_dimension pft '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'levgrnd', dimid), & 'get_sparse_dims','inq_dimid levgrnd '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=mylevgrnd), & 'get_sparse_dims','inquire_dimension levgrnd '//trim(fname)) if (mylevgrnd /= nlevgrnd) then write(string1,*)'Number of ground levels in restart file is',mylevgrnd write(string2,*)'Number of ground levels in history file is', nlevgrnd call error_handler(E_ERR,'get_sparse_dims',string1,source,revision,revdate,text2=string2) endif call nc_check(nf90_inq_dimid(ncid, 'levlak', dimid), & 'get_sparse_dims','inq_dimid levlak '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevlak), & 'get_sparse_dims','inquire_dimension levlak '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'levtot', dimid), & 'get_sparse_dims','inq_dimid levtot '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevtot), & 'get_sparse_dims','inquire_dimension levtot '//trim(fname)) call nc_check(nf90_inq_dimid(ncid, 'numrad', dimid), & 'get_sparse_dims','inq_dimid numrad '//trim(fname)) call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nnumrad), & 'get_sparse_dims','inquire_dimension numrad '//trim(fname)) ! levsno is presently required, but I can envision a domain/experiment that ! will not have snow levels. How this relates to variables dimensioned 'levtot' ! is unclear. For that reason, levsno is presently required. istatus = nf90_inq_dimid(ncid, 'levsno', dimid) if (istatus == nf90_noerr) then call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevsno), & 'get_sparse_dims','inquire_dimension levsno '//trim(fname)) endif ! levsno1, rtmlon, rtmlat are optional ... so it is not a fatal error if they are not present. istatus = nf90_inq_dimid(ncid, 'levsno1', dimid) if (istatus == nf90_noerr) then call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevsno1), & 'get_sparse_dims','inquire_dimension levsno1 '//trim(fname)) endif istatus = nf90_inq_dimid(ncid, 'rtmlon', dimid) if (istatus == nf90_noerr) then call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nrtmlon), & 'get_sparse_dims','inquire_dimension rtmlon '//trim(fname)) endif istatus = nf90_inq_dimid(ncid, 'rtmlat', dimid) if (istatus == nf90_noerr) then call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nrtmlat), & 'get_sparse_dims','inquire_dimension rtmlat '//trim(fname)) endif if (cstat == 'close') then call nc_check(nf90_close(ncid),'get_sparse_dims','close '//trim(fname) ) ncid = 0 endif ! Echo what we know if desired. if ((debug > 7) .and. do_output()) then write(logfileunit,*) write(logfileunit,*)'get_sparse_dims output follows:' write(logfileunit,*)'Ngridcell = ',Ngridcell write(logfileunit,*)'Nlandunit = ',Nlandunit write(logfileunit,*)'Ncolumn = ',Ncolumn write(logfileunit,*)'Npft = ',Npft write(logfileunit,*)'Nlevgrnd = ',Nlevgrnd write(logfileunit,*)'Nlevlak = ',Nlevlak write(logfileunit,*)'Nlevsno = ',Nlevsno write(logfileunit,*)'Nlevsno1 = ',Nlevsno1 write(logfileunit,*)'Nlevtot = ',Nlevtot write(logfileunit,*)'Nnumrad = ',Nnumrad write(logfileunit,*)'Nrtmlon = ',Nrtmlon write(logfileunit,*)'Nrtmlat = ',Nrtmlat write( * ,*) write( * ,*)'get_sparse_dims output follows:' write( * ,*)'Ngridcell = ',Ngridcell write( * ,*)'Nlandunit = ',Nlandunit write( * ,*)'Ncolumn = ',Ncolumn write( * ,*)'Npft = ',Npft write( * ,*)'Nlevgrnd = ',Nlevgrnd write( * ,*)'Nlevlak = ',Nlevlak write( * ,*)'Nlevsno = ',Nlevsno write( * ,*)'Nlevsno1 = ',Nlevsno1 write( * ,*)'Nlevtot = ',Nlevtot write( * ,*)'Nnumrad = ',Nnumrad write( * ,*)'Nrtmlon = ',Nrtmlon write( * ,*)'Nrtmlat = ',Nrtmlat endif end subroutine get_sparse_dims !------------------------------------------------------------------ subroutine get_sparse_geog(ncid, fname, cstat) !------------------------------------------------------------------ ! ! Read the geography information from from the restart netcdf file. integer, intent(inout) :: ncid character(len=*), intent(in) :: fname character(len=*), intent(in) :: cstat integer :: VarID if (ncid == 0) then call nc_check(nf90_open(trim(fname), nf90_nowrite, ncid), & 'get_sparse_geog','open '//trim(fname)) endif ! Make sure the variables are the right size ... ! by comparing agains the size of the variable ... if ( Ngridcell < 0 ) then write(string1,*)'Unable to read the number of gridcells.' call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate) endif if ( Nlandunit < 0 ) then write(string1,*)'Unable to read the number of land units.' call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate) endif if ( Ncolumn < 0 ) then write(string1,*)'Unable to read the number of columns.' call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate) endif if ( Npft < 0 ) then write(string1,*)'Unable to read the number of pfts.' call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate) endif ! Read the netcdf file data call nc_check(nf90_inq_varid(ncid, 'grid1d_ixy', VarID), 'get_sparse_geog', & 'inq_varid grid1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, grid1d_ixy), 'get_sparse_geog', & 'get_var grid1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'grid1d_jxy', VarID), 'get_sparse_geog', & 'inq_varid grid1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, grid1d_jxy), 'get_sparse_geog', & 'get_var grid1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'land1d_ixy', VarID), 'get_sparse_geog', & 'inq_varid land1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, land1d_ixy), 'get_sparse_geog', & 'get_var land1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'land1d_jxy', VarID), 'get_sparse_geog', & 'inq_varid land1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, land1d_jxy), 'get_sparse_geog', & 'get_var land1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'land1d_wtxy', VarID), 'get_sparse_geog', & 'inq_varid land1d_wtxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, land1d_wtxy), 'get_sparse_geog', & 'get_var land1d_wtxy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'cols1d_ixy', VarID), 'get_sparse_geog', & 'inq_varid cols1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, cols1d_ixy), 'get_sparse_geog', & 'get_var cols1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'cols1d_jxy', VarID), 'get_sparse_geog', & 'inq_varid cols1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, cols1d_jxy), 'get_sparse_geog', & 'get_var cols1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'cols1d_wtxy', VarID), 'get_sparse_geog', & 'inq_varid cols1d_wtxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, cols1d_wtxy), 'get_sparse_geog', & 'get_var cols1d_wtxy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'cols1d_ityplun', VarID), 'get_sparse_geog', & 'inq_varid cols1d_ityplun '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, cols1d_ityplun), 'get_sparse_geog', & 'get_var cols1d_ityplun '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'pfts1d_ixy', VarID), 'get_sparse_geog', & 'inq_varid pfts1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, pfts1d_ixy), 'get_sparse_geog', & 'get_var pfts1d_ixy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'pfts1d_jxy', VarID), 'get_sparse_geog', & 'inq_varid pfts1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, pfts1d_jxy), 'get_sparse_geog', & 'get_var pfts1d_jxy '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'pfts1d_wtxy', VarID), 'get_sparse_geog', & 'inq_varid pfts1d_wtxy '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, pfts1d_wtxy), 'get_sparse_geog', & 'get_var pfts1d_wtxy '//trim(clm_restart_filename)) ! zsno is NOT optional ... so it IS a fatal error if it is not present (for now, anyway). ! as read into fortran ... zsno(:,1) is the level closest to the sun. ! as read into fortran ... zsno(:,5) is the level closest to the ground. if (nlevsno > 0 ) then call nc_check(nf90_inq_varid(ncid, 'ZSNO', VarID), & & 'get_sparse_geog', 'inq_varid ZSNO '//trim(fname)) call nc_check(nf90_get_var( ncid, VarID, zsno), & & 'get_sparse_geog', 'get_var ZSNO '//trim(fname)) else write(string1,*) 'levsno must be in restart file' call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate) endif if (cstat == 'close') then call nc_check(nf90_close(ncid),'get_sparse_geog','close '//trim(fname) ) ncid = 0 endif ! A little sanity check if ((debug > 7) .and. do_output()) then write(logfileunit,*) write(logfileunit,*)'Raw lat/lon information as read ...' write(logfileunit,*)'grid1d_ixy range ',minval(grid1d_ixy), maxval(grid1d_ixy) write(logfileunit,*)'grid1d_jxy range ',minval(grid1d_jxy), maxval(grid1d_jxy) write(logfileunit,*)'land1d_ixy range ',minval(land1d_ixy), maxval(land1d_ixy) write(logfileunit,*)'land1d_jxy range ',minval(land1d_jxy), maxval(land1d_jxy) write(logfileunit,*)'land1d_wtxy range ',minval(land1d_wtxy), maxval(land1d_wtxy) write(logfileunit,*)'cols1d_ixy range ',minval(cols1d_ixy), maxval(cols1d_ixy) write(logfileunit,*)'cols1d_jxy range ',minval(cols1d_jxy), maxval(cols1d_jxy) write(logfileunit,*)'cols1d_wtxy range ',minval(cols1d_wtxy), maxval(cols1d_wtxy) write(logfileunit,*)'cols1d_ityplun range ',minval(cols1d_ityplun),maxval(cols1d_ityplun) write(logfileunit,*)'pfts1d_ixy range ',minval(pfts1d_ixy), maxval(pfts1d_ixy) write(logfileunit,*)'pfts1d_jxy range ',minval(pfts1d_jxy), maxval(pfts1d_jxy) write(logfileunit,*)'pfts1d_wtxy range ',minval(pfts1d_wtxy), maxval(pfts1d_wtxy) if (nlevsno > 0) write(logfileunit,*)'zsno range ',minval(zsno),maxval(zsno) write( * ,*) write( * ,*)'Raw lat/lon information as read ...' write( * ,*)'grid1d_ixy range ',minval(grid1d_ixy), maxval(grid1d_ixy) write( * ,*)'grid1d_jxy range ',minval(grid1d_jxy), maxval(grid1d_jxy) write( * ,*)'land1d_ixy range ',minval(land1d_ixy), maxval(land1d_ixy) write( * ,*)'land1d_jxy range ',minval(land1d_jxy), maxval(land1d_jxy) write( * ,*)'land1d_wtxy range ',minval(land1d_wtxy), maxval(land1d_wtxy) write( * ,*)'cols1d_ixy range ',minval(cols1d_ixy), maxval(cols1d_ixy) write( * ,*)'cols1d_jxy range ',minval(cols1d_jxy), maxval(cols1d_jxy) write( * ,*)'cols1d_wtxy range ',minval(cols1d_wtxy), maxval(cols1d_wtxy) write( * ,*)'cols1d_ityplun range ',minval(cols1d_ityplun),maxval(cols1d_ityplun) write( * ,*)'pfts1d_ixy range ',minval(pfts1d_ixy), maxval(pfts1d_ixy) write( * ,*)'pfts1d_jxy range ',minval(pfts1d_jxy), maxval(pfts1d_jxy) write( * ,*)'pfts1d_wtxy range ',minval(pfts1d_wtxy), maxval(pfts1d_wtxy) if (nlevsno > 0) write( * ,*)'zsno range ',minval(zsno),maxval(zsno) endif return end subroutine get_sparse_geog !------------------------------------------------------------------ function get_state_time_ncid( ncid ) !------------------------------------------------------------------ ! The restart netcdf files have the time of the state. type(time_type) :: get_state_time_ncid integer, intent(in) :: ncid integer :: VarID integer :: rst_curr_ymd, rst_curr_tod, leftover integer :: year, month, day, hour, minute, second if ( .not. module_initialized ) call static_init_model call nc_check(nf90_inq_varid(ncid, 'timemgr_rst_curr_ymd', VarID), 'get_state_time_ncid', & & 'inq_varid timemgr_rst_curr_ymd '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, rst_curr_ymd), 'get_state_time_ncid', & & 'get_var rst_curr_ymd '//trim(clm_restart_filename)) call nc_check(nf90_inq_varid(ncid, 'timemgr_rst_curr_tod', VarID), 'get_state_time_ncid', & & 'inq_varid timemgr_rst_curr_tod '//trim(clm_restart_filename)) call nc_check(nf90_get_var( ncid, VarID, rst_curr_tod), 'get_state_time_ncid', & & 'get_var rst_curr_tod '//trim(clm_restart_filename)) year = rst_curr_ymd/10000 leftover = rst_curr_ymd - year*10000 month = leftover/100 day = leftover - month*100 hour = rst_curr_tod/3600 leftover = rst_curr_tod - hour*3600 minute = leftover/60 second = leftover - minute*60 get_state_time_ncid = set_date(year, month, day, hour, minute, second) end function get_state_time_ncid !------------------------------------------------------------------ function get_state_time_fname(filename) !------------------------------------------------------------------ ! the static_init_model ensures that the clm namelists are read. ! type(time_type) :: get_state_time_fname character(len=*), intent(in) :: filename integer :: ncid if ( .not. module_initialized ) call static_init_model if ( .not. file_exist(filename) ) then write(string1,*) 'cannot open file ', trim(filename),' for reading.' call error_handler(E_ERR,'get_state_time_fname',string1,source,revision,revdate) endif call nc_check( nf90_open(trim(filename), NF90_NOWRITE, ncid), & 'get_state_time_fname', 'open '//trim(filename)) get_state_time_fname = get_state_time_ncid(ncid) call nc_check(nf90_close(ncid),'get_state_time_fname', 'close '//trim(filename)) end function get_state_time_fname !------------------------------------------------------------------ function set_model_time_step() !------------------------------------------------------------------ ! This defines the window used for assimilation. ! all observations +/- half this timestep are assimilated. type(time_type) :: set_model_time_step set_model_time_step = set_time(assimilation_period_seconds, assimilation_period_days) end function set_model_time_step !------------------------------------------------------------------ subroutine get_clm_restart_filename( filename ) character(len=*), intent(OUT) :: filename if ( .not. module_initialized ) call static_init_model filename = trim(clm_restart_filename) end subroutine get_clm_restart_filename subroutine get_clm_history_filename( filename ) character(len=*), intent(OUT) :: filename if ( .not. module_initialized ) call static_init_model filename = trim(clm_history_filename) end subroutine get_clm_history_filename !------------------------------------------------------------------ !> This routine checks the user input against the variables available in the !> input netcdf file to see if it is possible to construct the DART state vector !> specified by the input.nml:model_nml:clm_variables variable. !> Each variable must have 6 entries. !> 1: variable name !> 2: DART KIND !> 3: minimum value - as a character string - if none, use 'NA' !> 4: maximum value - as a character string - if none, use 'NA' !> 5: what file contains the variable - '.r. => restart', '.h0. => h0 history file' !> 6: does the variable get updated in the restart file or not ... !> only variables from restart files may be updated. !> 'UPDATE' => update the variable in the restart file !> 'NO_COPY_BACK' => do not copy the variable back to the restart file !> all these variables will be updated INTERNALLY IN DART !> only variables marked '.r', 'UPDATE' will be modified for CLM. !> !> The calling code should check to see if the variable exists. subroutine parse_variable_table( state_variables, ngood, table ) character(len=*), dimension(:), intent(in) :: state_variables integer, intent(out) :: ngood character(len=*), dimension(:,:), intent(out) :: table integer :: nrows, ncols, i character(len=NF90_MAX_NAME) :: varname ! column 1 character(len=NF90_MAX_NAME) :: dartstr ! column 2 character(len=NF90_MAX_NAME) :: minvalstring ! column 3 character(len=NF90_MAX_NAME) :: maxvalstring ! column 4 character(len=NF90_MAX_NAME) :: origin_file ! column 5 character(len=NF90_MAX_NAME) :: state_or_aux ! column 6 nrows = size(table,1) ncols = size(table,2) ! This loop just repackages the 1D array of values into a 2D array. ! We can do some miniminal checking along the way. ! Determining which file to check is going to be more complicated. ngood = 0 MyLoop : do i = 1, nrows varname = trim(state_variables(ncols*i - 5)) dartstr = trim(state_variables(ncols*i - 4)) minvalstring = trim(state_variables(ncols*i - 3)) maxvalstring = trim(state_variables(ncols*i - 2)) origin_file = trim(state_variables(ncols*i - 1)) state_or_aux = trim(state_variables(ncols*i )) call to_upper(origin_file) call to_upper(state_or_aux) table(i,1) = trim(varname) table(i,2) = trim(dartstr) table(i,3) = trim(minvalstring) table(i,4) = trim(maxvalstring) table(i,5) = trim(origin_file) table(i,6) = trim(state_or_aux) ! If the first element is empty, we have found the end of the list. if ( table(i,1) == ' ' ) exit MyLoop ! Any other condition is an error. if ( any(table(i,:) == ' ') ) then string1 = 'input.nml &model_nml:clm_variables not fully specified' string2 = 'must be 6 entries per variable. Last known variable name is' string3 = '['//trim(table(i,1))//'] ... (without the [], naturally)' call error_handler(E_ERR, 'parse_variable_table', string1, & source, revision, revdate, text2=string2, text3=string3) endif ! Make sure DART kind is valid if( get_raw_obs_kind_index(dartstr) < 0 ) then write(string1,'(''there is no obs_kind <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) call error_handler(E_ERR,'parse_variable_table',string1,source,revision,revdate) endif ! Record the contents of the DART state vector if ((debug > 8) .and. do_output()) then write(logfileunit,*)'variable ',i,' is ',trim(table(i,1)), ' ', trim(table(i,2)),' ', & trim(table(i,3)), ' ', trim(table(i,4)),' ', & trim(table(i,5)), ' ', trim(table(i,6)) write( * ,*)'variable ',i,' is ',trim(table(i,1)), ' ', trim(table(i,2)),' ', & trim(table(i,3)), ' ', trim(table(i,4)),' ', & trim(table(i,5)), ' ', trim(table(i,6)) endif ngood = ngood + 1 enddo MyLoop if (ngood == nrows) then string1 = 'WARNING: There is a possibility you need to increase ''max_state_variables''' write(string2,'(''WARNING: you have specified at least '',i4,'' perhaps more.'')')ngood call error_handler(E_MSG,'parse_variable_table',string1,text2=string2) endif ! Check to see if zsno is part of the requested variables end subroutine parse_variable_table !------------------------------------------------------------------ function FindTimeDimension(ncid) result(timedimid) ! Find the Time Dimension ID in a netCDF file. ! If there is none - (spelled the obvious way) - the routine ! returns a negative number. You don't HAVE to have a TIME dimension. integer :: timedimid integer, intent(in) :: ncid integer :: nc_rc TimeDimID = -1 ! same as the netCDF library routines. nc_rc = nf90_inq_dimid(ncid,'TIME',dimid=TimeDimID) if ( nc_rc /= NF90_NOERR ) then ! did not find it - try another spelling nc_rc = nf90_inq_dimid(ncid,'Time',dimid=TimeDimID) if ( nc_rc /= NF90_NOERR ) then ! did not find it - try another spelling nc_rc = nf90_inq_dimid(ncid,'time',dimid=TimeDimID) endif endif end function FindTimeDimension !------------------------------------------------------------------ !> define_var_dims() takes the N-dimensional variable and appends the DART !> dimensions of 'copy' and 'time'. If the variable initially had a 'time' !> dimension, it is ignored because (by construction) it is a singleton !> dimension. subroutine define_var_dims(ivar, ncid, memberdimid, unlimiteddimid, ndims, dimids) integer, intent(in) :: ivar, ncid, memberdimid, unlimiteddimid integer, intent(out) :: ndims integer, dimension(:), intent(out) :: dimids character(len=NF90_MAX_NAME),dimension(NF90_MAX_VAR_DIMS) :: dimnames integer :: i, mydimid ndims = 0 DIMLOOP : do i = 1,progvar(ivar)%numdims if (progvar(ivar)%dimnames(i) == 'time') cycle DIMLOOP call nc_check(nf90_inq_dimid(ncid=ncid, name=progvar(ivar)%dimnames(i), dimid=mydimid), & 'define_var_dims','inq_dimid '//trim(progvar(ivar)%dimnames(i))) ndims = ndims + 1 dimids(ndims) = mydimid dimnames(ndims) = progvar(ivar)%dimnames(i) enddo DIMLOOP ! The last two dimensions are always 'copy' and 'time' ndims = ndims + 1 dimids(ndims) = memberdimid dimnames(ndims) = 'copy' ndims = ndims + 1 dimids(ndims) = unlimitedDimid dimnames(ndims) = 'time' if ((debug > 8) .and. do_output()) then write(logfileunit,*) write(logfileunit,*)'define_var_dims knowledge' write(logfileunit,*)trim(progvar(ivar)%varname),' has original dimnames: ', & (/( trim(progvar(ivar)%dimnames(i))//' ',i=1,progvar(ivar)%numdims) /) write(logfileunit,*)trim(progvar(ivar)%varname),' repackaging into dimnames: ', & (/ (trim(dimnames(i))//' ',i=1,ndims) /) write(logfileunit,*)'thus dimids ',dimids(1:ndims) write( * ,*) write( * ,*)'define_var_dims knowledge' write( * ,*)trim(progvar(ivar)%varname),' has original dimnames: ', & (/( trim(progvar(ivar)%dimnames(i))//' ',i=1,progvar(ivar)%numdims) /) write( * ,*)trim(progvar(ivar)%varname),' repackaging into dimnames: ', & (/ (trim(dimnames(i))//' ',i=1,ndims) /) write( * ,*)'thus dimids ',dimids(1:ndims) endif return end subroutine define_var_dims !------------------------------------------------------------------ subroutine fill_levels(dimname,icol,enlevels,levtot) ! subroutine fill_levels(dimname,icol,enlevels,levtot) ! ! dimname ... is it dimensioned 'levgrnd' or 'levsno' or 'levtot' ... ! icol ... which CLM 'column' are we in ! enlevels ... the expected number of levels ... varshape ! levtot ... the output array of vertical coordinates ! ! The total number of levels is defined to be the soil levels (fixed) ! plus the number of snow levels, which can vary by column. ! The history file contains the depths of the soil levels (ncdf var 'levgrnd'); ! these are in levtot(1:nlevgrnd). ! The restart file contains the depths of the snow levels (ncdf var 'ZSNO'). ! The tricky bit is that they are in reverse order ... and differ from one model to another. ! If you simply grab the netcdf variable zsno, ! the level closest to the soil is actually the highest index. ! ! From Matlab (which indexes like Fortran) !> size(zsno) ans = 13693 5 !> zsno(1,:) ans = -1.4202 -1.3852 -1.3052 -1.1352 -0.5101 ! | | | | |...... closest to soil surface ! | | | |............... one level 'up' ! | | |......................... one level 'up' ! | |................................... one level 'up' ! |.............................................. closest to sun ! ! If there is no snow ... the corresponding zsno is ZERO ... !> zsno(508,:) ans = 0 -0.5736 -0.5389 -0.4591 -0.2021 ! ! The following Matlab code may be used to explore a variable's storage order. ! (a better version is in the clm/matlab/CheckStorageOrder.m function) ! ! h2o = nc_varget(fname,'H2OSOI_LIQ'); ! h2o(h2o > 1.0E30) = NaN; ! lat = nc_varget(fname,'cols1d_lat'); ! lon = nc_varget(fname,'cols1d_lon'); ! figure(1); plot3(lon,lat,h2o(:,1),'x'); hold on; worldmap; view(0,90) ! figure(2); plot3(lon,lat,h2o(:,2),'x'); hold on; worldmap; view(0,90) ! figure(3); plot3(lon,lat,h2o(:,3),'x'); hold on; worldmap; view(0,90) ! figure(4); plot3(lon,lat,h2o(:,4),'x'); hold on; worldmap; view(0,90) ! figure(5); plot3(lon,lat,h2o(:,5),'x'); hold on; worldmap; view(0,90) ! figure(6); plot3(lon,lat,h2o(:,6),'x'); hold on; worldmap; view(0,90) character(len=*), intent(in) :: dimname integer, intent(in) :: icol integer, intent(in) :: enlevels real(r8), dimension(:), intent(inout) :: levtot if (dimname == 'levsno') then if (nlevsno /= enlevels) then write(string1,*) 'dimension ', trim(dimname),' has declared length ',enlevels write(string2,*) 'not the known number of snow levels ',nlevsno call error_handler(E_ERR,'fill_levels', string1, & source, revision, revdate, text2=string2) endif levtot(1:nlevsno) = zsno(1:nlevsno,icol) elseif (dimname == 'levgrnd') then if (nlevgrnd /= enlevels) then write(string1,*) 'dimension ', trim(dimname),' has declared length ',enlevels write(string2,*) 'not the known number of soil levels ',nlevgrnd call error_handler(E_ERR,'fill_levels', string1, & source, revision, revdate, text2=string2) endif levtot(1:nlevgrnd) = LEVGRND elseif (dimname == 'levtot') then ! This block assumes anything dimensioned 'levtot' has the first nlevsno levels ! followed by nlevgrnd levels. Dunno what to do with lake stuff ... if (nlevtot /= enlevels) then write(string1,*) 'dimension ', trim(dimname),' has declared length ',enlevels write(string2,*) 'not the known number of total levels ',nlevtot call error_handler(E_ERR,'fill_levels', string1, & source, revision, revdate, text2=string2) endif if (nlevtot /= nlevgrnd + nlevsno) then write(string1,*) 'nlevtot ', nlevtot,' is not equal to nlevgrnd + nlevsno' write(string2,*) 'nlevgrnd is ',nlevgrnd,' nlevsno is ',nlevsno,' total of ',nlevgrnd+nlevsno call error_handler(E_ERR,'fill_levels', string1, & source, revision, revdate, text2=string2) endif levtot(1:nlevsno) = zsno(1:nlevsno,icol) levtot(nlevsno+1:nlevsno+nlevgrnd) = LEVGRND else write(string1,*) 'Unable to determine vertical coordinate for column ',icol write(string2,*) 'unknown dimension name: ',trim(dimname) call error_handler(E_ERR,'fill_levels', string1, & source, revision, revdate, text2=string2) endif end subroutine fill_levels !------------------------------------------------------------------ subroutine gridcell_components(varstring) ! In order to exercise some of the routines, it is necessary to know ! which gridcells have multiple land units ! which gridcells have multiple columns ! which gridcells have multiple PFTs ! ! This routine simply tells me which gridcells are 'interesting'. ! Each level counts separately. 1 column with 20 levels ... yields a count of 20 ! ! It is very similar to SetLocatorArrays(), but only does the landunit,column,pft ! that the variable uses. It is also public - currently only used by ! model_mod_check. character(len=*), intent(in) :: varstring ! T_SOISNO, H2OSOI_LIQ ! Local storage integer :: ivar, indexi, i, j integer, allocatable, dimension(:,:) :: countmat if ( .not. module_initialized ) call static_init_model allocate(countmat(nlon,nlat)) countmat = 0 VARTYPES : do ivar = 1,nfields ! Skip to the right variable if ( trim(progvar(ivar)%varname) /= varstring) cycle VARTYPES ! Create a count of all the multiples in a gridcell do indexi = progvar(ivar)%index1, progvar(ivar)%indexN i = lonixy(indexi) j = latjxy(indexi) countmat(i,j) = countmat(i,j) + 1 enddo enddo VARTYPES if (do_output() .and. (debug > 0)) then write(*,*)'exploring '//trim(varstring) do j = 1,nlat do i = 1,nlon if ( countmat(i,j) > 1) & write(*,'(''gridcell'',2(1x,i8),'' has '',i6,'' lon/lat'',2(1x,f12.7))') & i,j,countmat(i,j),LON(i),LAT(j) enddo enddo endif deallocate(countmat) end subroutine gridcell_components subroutine get_var_1d(ncid, varname, var1d) ! This function will return a R8 array with the netCDF attributes applied. ! scale_factor, offset will be applied, ! missing_value, _FillValue will be replaced by the DART missing value ... ! If _FillValue is defined then it should be scalar and of the same type as the variable. ! If the variable is packed using scale_factor and add_offset attributes (see below), ! the _FillValue attribute should have the data type of the packed data. ! ! missing_value ! When scale_factor and add_offset are used for packing, the value(s) of the missing_value ! attribute should be specified in the domain of the data in the file (the packed data), ! so that missing values can be detected before the scale_factor and add_offset are applied. ! ! scale_factor ! If present for a variable, the data are to be multiplied by this factor after the data ! are read by the application that accesses the data. If valid values are specified using ! the valid_min, valid_max, valid_range, or _FillValue attributes, those values should be ! specified in the domain of the data in the file (the packed data), so that they can be ! interpreted before the scale_factor and add_offset are applied. ! ! add_offset ! If present for a variable, this number is to be added to the data after it is read by ! the application that accesses the data. If both scale_factor and add_offset attributes ! are present, the data are first scaled before the offset is added. integer, intent(in) :: ncid character(len=*), intent(in) :: varname real(r8), dimension(:), intent(out) :: var1d integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, dimlens, ncstart, nccount integer :: VarID, numdims, xtype, io1, io2 integer :: TimeDimID, time_dimlen, timeindex integer :: spvalINT real(r4) :: spvalR4 real(r8) :: scale_factor, add_offset real(r8) :: spvalR8 integer, allocatable, dimension(:) :: intarray real(r4), allocatable, dimension(:) :: r4array if ( .not. module_initialized ) call static_init_model ! a little whitespace makes this a lot more readable if (do_output() .and. (debug > 1)) then write(*,*) write(logfileunit,*) endif io1 = nf90_inq_dimid(ncid, 'time', TimeDimID) if (io1 /= NF90_NOERR) TimeDimID = MISSING_I call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), 'get_var_1d', 'inq_varid '//varname) call nc_check(nf90_inquire_variable( ncid, VarID, dimids=dimIDs, ndims=numdims, & xtype=xtype), 'get_var_1d', 'inquire_variable '//varname) call nc_check(nf90_inquire_dimension(ncid, dimIDs(1), len=dimlens(1)), & 'get_var_1d', 'inquire_dimension '//varname) if ((numdims /= 1) .or. (size(var1d) /= dimlens(1)) ) then write(string1,*) trim(varname)//' is not the expected shape/length of ', size(var1d) call error_handler(E_ERR,'get_var_1d',string1,source,revision,revdate) endif ncstart = 1 nccount = dimlens(1) if (dimIDs(1) == TimeDimID) then call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=time_dimlen ), & 'get_var_1d', 'inquire_dimension time '//trim(varname)) timeindex = FindDesiredTimeIndx(ncid, time_dimlen, varname) ncstart(1) = timeindex nccount(1) = 1 endif if (do_output() .and. (debug > 1)) then write(*,*)'get_var_1d: variable ['//trim(varname)//']' write(*,*)'get_var_1d: start ',ncstart(1:numdims) write(*,*)'get_var_1d: count ',nccount(1:numdims) endif if (xtype == NF90_INT) then allocate(intarray(dimlens(1))) call nc_check(nf90_get_var(ncid, VarID, values=intarray, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_1d', 'get_var '//varname) var1d = intarray ! perform type conversion to desired output type io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalINT) if ( io1 == NF90_NOERR) where (intarray == spvalINT) var1d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalINT,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_1d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalINT) if ( io2 == NF90_NOERR) where (intarray == spvalINT) var1d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalINT,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_1d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var1d /= MISSING_R8) var1d = var1d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var1d /= MISSING_R8) var1d = var1d * scale_factor elseif (io2 == NF90_NOERR) then where (var1d /= MISSING_R8) var1d = var1d + add_offset endif deallocate(intarray) elseif (xtype == NF90_FLOAT) then allocate(r4array(dimlens(1))) call nc_check(nf90_get_var(ncid, VarID, values=r4array, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_1d', 'get_var '//varname) var1d = r4array ! perform type conversion to desired output type io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalR4) if ( io1 == NF90_NOERR) where (r4array == spvalR4) var1d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalR4,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_1d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalR4) if ( io2 == NF90_NOERR) where (r4array == spvalR4) var1d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalR4,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_1d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var1d /= MISSING_R8) var1d = var1d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var1d /= MISSING_R8) var1d = var1d * scale_factor elseif (io2 == NF90_NOERR) then where (var1d /= MISSING_R8) var1d = var1d + add_offset endif deallocate(r4array) elseif (xtype == NF90_DOUBLE) then call nc_check(nf90_get_var(ncid, VarID, values=var1d, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_1d', 'get_var '//varname) io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalR8) if ( io1 == NF90_NOERR) where (var1d == spvalR8) var1d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalR8,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_1d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalR8) if ( io2 == NF90_NOERR) where (var1d == spvalR8) var1d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalR8,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_1d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var1d /= MISSING_R8) var1d = var1d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var1d /= MISSING_R8) var1d = var1d * scale_factor elseif (io2 == NF90_NOERR) then where (var1d /= MISSING_R8) var1d = var1d + add_offset endif else write(string1,*) trim(varname)//' has unsupported (by DART) xtype of', xtype call error_handler(E_ERR,'get_var_1d',string1,source,revision,revdate) endif end subroutine get_var_1d subroutine get_var_2d(ncid, varname, var2d) ! This function will return a R8 array with the netCDF attributes applied. ! scale_factor, offset will be applied, ! missing_value, _FillValue will be replaced by the DART missing value ... ! If _FillValue is defined then it should be scalar and of the same type as the variable. ! If the variable is packed using scale_factor and add_offset attributes (see below), ! the _FillValue attribute should have the data type of the packed data. ! ! missing_value ! When scale_factor and add_offset are used for packing, the value(s) of the missing_value ! attribute should be specified in the domain of the data in the file (the packed data), ! so that missing values can be detected before the scale_factor and add_offset are applied. ! ! scale_factor ! If present for a variable, the data are to be multiplied by this factor after the data ! are read by the application that accesses the data. If valid values are specified using ! the valid_min, valid_max, valid_range, or _FillValue attributes, those values should be ! specified in the domain of the data in the file (the packed data), so that they can be ! interpreted before the scale_factor and add_offset are applied. ! ! add_offset ! If present for a variable, this number is to be added to the data after it is read by ! the application that accesses the data. If both scale_factor and add_offset attributes ! are present, the data are first scaled before the offset is added. integer, intent(in) :: ncid character(len=*), intent(in) :: varname real(r8), dimension(:,:), intent(out) :: var2d integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, dimlens, ncstart, nccount integer :: VarID, numdims, xtype, io1, io2, i integer :: TimeDimID, time_dimlen, timeindex integer :: spvalINT real(r4) :: spvalR4 real(r8) :: scale_factor, add_offset real(r8) :: spvalR8 integer, allocatable, dimension(:,:) :: intarray real(r4), allocatable, dimension(:,:) :: r4array if ( .not. module_initialized ) call static_init_model ! a little whitespace makes this a lot more readable if (do_output() .and. (debug > 1)) then write(*,*) write(logfileunit,*) endif io1 = nf90_inq_dimid(ncid, 'time', TimeDimID) if (io1 /= NF90_NOERR) TimeDimID = MISSING_I call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), 'get_var_2d', 'inq_varid') call nc_check(nf90_inquire_variable( ncid, VarID, dimids=dimIDs, ndims=numdims, & xtype=xtype), 'get_var_2d', 'inquire_variable') if ( (numdims /= 2) ) then write(string1,*) trim(varname)//' is not a 2D variable as expected.' call error_handler(E_ERR,'get_var_2d',string1,source,revision,revdate) endif ncstart(:) = 1 nccount(:) = 1 DimCheck : do i = 1,numdims write(string1,'(a,i2,1x,A)') 'inquire dimension ',i,trim(varname) call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlens(i)), & 'get_var_2d', string1) if ( dimIDs(i) == TimeDimID ) then call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=time_dimlen), & 'get_var_2d', 'inquire_dimension time '//trim(varname)) timeindex = FindDesiredTimeIndx(ncid, time_dimlen, varname) ncstart(i) = timeindex dimlens(i) = 1 elseif ( size(var2d,i) /= dimlens(i) ) then write(string1,*) trim(varname)//' has shape ', dimlens(1:numdims) write(string2,*) 'which is not the expected shape of ', size(var2d,1),size(var2d,2) call error_handler(E_ERR,'get_var_2d',string1,source,revision,revdate,text2=string2) endif nccount(i) = dimlens(i) enddo DimCheck if (do_output() .and. (debug > 1)) then write(*,*)'get_var_2d: variable ['//trim(varname)//']' write(*,*)'get_var_2d: start ',ncstart(1:numdims) write(*,*)'get_var_2d: count ',nccount(1:numdims) endif if (xtype == NF90_INT) then allocate(intarray(dimlens(1),dimlens(2))) call nc_check(nf90_get_var(ncid, VarID, values=intarray, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_2d', 'get_var '//varname) var2d = intarray ! perform type conversion to desired output type io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalINT) if ( io1 == NF90_NOERR) where (intarray == spvalINT) var2d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalINT,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_2d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalINT) if ( io2 == NF90_NOERR) where (intarray == spvalINT) var2d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalINT,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_2d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var2d /= MISSING_R8) var2d = var2d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var2d /= MISSING_R8) var2d = var2d * scale_factor elseif (io2 == NF90_NOERR) then where (var2d /= MISSING_R8) var2d = var2d + add_offset endif deallocate(intarray) elseif (xtype == NF90_FLOAT) then allocate(r4array(dimlens(1),dimlens(2))) call nc_check(nf90_get_var(ncid, VarID, values=r4array, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_2d', 'get_var '//varname) var2d = r4array ! perform type conversion to desired output type io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalR4) if ( io1 == NF90_NOERR) where (r4array == spvalR4) var2d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalR4,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_2d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalR4) if ( io2 == NF90_NOERR) where (r4array == spvalR4) var2d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalR4,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_2d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var2d /= MISSING_R8) var2d = var2d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var2d /= MISSING_R8) var2d = var2d * scale_factor elseif (io2 == NF90_NOERR) then where (var2d /= MISSING_R8) var2d = var2d + add_offset endif deallocate(r4array) elseif (xtype == NF90_DOUBLE) then call nc_check(nf90_get_var(ncid, VarID, values=var2d, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_2d', 'get_var '//varname) io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalR8) if ( io1 == NF90_NOERR) where (var2d == spvalR8) var2d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalR8,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_2d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalR8) if ( io2 == NF90_NOERR) where (var2d == spvalR8) var2d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalR8,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_2d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var2d /= MISSING_R8) var2d = var2d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var2d /= MISSING_R8) var2d = var2d * scale_factor elseif (io2 == NF90_NOERR) then where (var2d /= MISSING_R8) var2d = var2d + add_offset endif else write(string1,*) trim(varname)//' has unsupported (by DART) xtype of', xtype call error_handler(E_ERR,'get_var_2d',string1,source,revision,revdate) endif end subroutine get_var_2d subroutine get_var_3d(ncid, varname, var3d) ! This function will return a R8 array with the netCDF attributes applied. ! scale_factor, offset will be applied, ! missing_value, _FillValue will be replaced by the DART missing value ... ! If _FillValue is defined then it should be scalar and of the same type as the variable. ! If the variable is packed using scale_factor and add_offset attributes (see below), ! the _FillValue attribute should have the data type of the packed data. ! ! missing_value ! When scale_factor and add_offset are used for packing, the value(s) of the missing_value ! attribute should be specified in the domain of the data in the file (the packed data), ! so that missing values can be detected before the scale_factor and add_offset are applied. ! ! scale_factor ! If present for a variable, the data are to be multiplied by this factor after the data ! are read by the application that accesses the data. If valid values are specified using ! the valid_min, valid_max, valid_range, or _FillValue attributes, those values should be ! specified in the domain of the data in the file (the packed data), so that they can be ! interpreted before the scale_factor and add_offset are applied. ! ! add_offset ! If present for a variable, this number is to be added to the data after it is read by ! the application that accesses the data. If both scale_factor and add_offset attributes ! are present, the data are first scaled before the offset is added. integer, intent(in) :: ncid character(len=*), intent(in) :: varname real(r8), dimension(:,:,:), intent(out) :: var3d integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, dimlens, ncstart, nccount integer :: i, TimeDimID, time_dimlen, timeindex integer :: VarID, numdims, xtype, io1, io2 integer :: spvalINT real(r4) :: spvalR4 real(r8) :: scale_factor, add_offset real(r8) :: spvalR8 integer, allocatable, dimension(:,:,:) :: intarray real(r4), allocatable, dimension(:,:,:) :: r4array if ( .not. module_initialized ) call static_init_model ! a little whitespace makes this a lot more readable if (do_output() .and. (debug > 1)) then write(*,*) write(logfileunit,*) endif ! 3D fields must have a time dimension. ! Need to know the Time Dimension ID and length call nc_check(nf90_inq_dimid(ncid, 'time', TimeDimID), & 'get_var_3d', 'inq_dimid time '//varname) call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=time_dimlen ), & 'get_var_3d', 'inquire_dimension time '//varname) call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), 'get_var_3d', 'inq_varid') call nc_check(nf90_inquire_variable( ncid, VarID, dimids=dimIDs, ndims=numdims, & xtype=xtype), 'get_var_3d', 'inquire_variable '//varname) if ( (numdims /= 3) ) then write(string1,*) trim(varname)//' is not a 3D variable as expected.' call error_handler(E_ERR,'get_var_3d',string1,source,revision,revdate) endif ! only expecting [nlon,nlat,time] ncstart(:) = 1 nccount(:) = 1 DimCheck : do i = 1,numdims write(string1,'(a,i2,1x,A)') 'inquire dimension ',i,trim(varname) call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlens(i)), & 'get_var_3d', trim(string1)) if ( dimIDs(i) == TimeDimID ) then timeindex = FindDesiredTimeIndx(ncid, time_dimlen, varname) ncstart(i) = timeindex dimlens(i) = 1 elseif ( size(var3d,i) /= dimlens(i) ) then write(string1,*) trim(varname)//' has shape ', dimlens(1:numdims) write(string2,*) 'which is not the expected shape of ', size(var3d,i) call error_handler(E_ERR,'get_var_3d',string1,source,revision,revdate,text2=string2) endif nccount(i) = dimlens(i) enddo DimCheck if (do_output() .and. (debug > 1)) then write(*,*)'get_var_3d: variable ['//trim(varname)//']' write(*,*)'get_var_3d: start ',ncstart(1:numdims) write(*,*)'get_var_3d: count ',nccount(1:numdims) endif if (xtype == NF90_INT) then allocate(intarray(dimlens(1),dimlens(2),dimlens(3))) call nc_check(nf90_get_var(ncid, VarID, values=intarray, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_3d', 'get_var '//varname) var3d = intarray ! perform type conversion to desired output type io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalINT) if ( io1 == NF90_NOERR) where (intarray == spvalINT) var3d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalINT,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_3d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalINT) if ( io2 == NF90_NOERR) where (intarray == spvalINT) var3d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalINT,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_3d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var3d /= MISSING_R8) var3d = var3d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var3d /= MISSING_R8) var3d = var3d * scale_factor elseif (io2 == NF90_NOERR) then where (var3d /= MISSING_R8) var3d = var3d + add_offset endif deallocate(intarray) elseif (xtype == NF90_FLOAT) then allocate(r4array(dimlens(1),dimlens(2),dimlens(3))) call nc_check(nf90_get_var(ncid, VarID, values=r4array, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_3d', 'get_var '//varname) var3d = r4array ! perform type conversion to desired output type io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalR4) if ( io1 == NF90_NOERR) where (r4array == spvalR4) var3d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalR4,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_3d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalR4) if ( io2 == NF90_NOERR) where (r4array == spvalR4) var3d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalR4,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_3d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var3d /= MISSING_R8) var3d = var3d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var3d /= MISSING_R8) var3d = var3d * scale_factor elseif (io2 == NF90_NOERR) then where (var3d /= MISSING_R8) var3d = var3d + add_offset endif deallocate(r4array) elseif (xtype == NF90_DOUBLE) then call nc_check(nf90_get_var(ncid, VarID, values=var3d, & start=ncstart(1:numdims), count=nccount(1:numdims)), & 'get_var_3d', 'get_var '//varname) io1 = nf90_get_att(ncid, VarID, '_FillValue' , spvalR8) if ( io1 == NF90_NOERR) where (var3d == spvalR8) var3d = MISSING_R8 if ( (io1 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing _FillValue ',spvalR8,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_3d',string1) endif io2 = nf90_get_att(ncid, VarID, 'missing_value' , spvalR8) if ( io2 == NF90_NOERR) where (var3d == spvalR8) var3d = MISSING_R8 if ( (io2 == NF90_NOERR) .and. do_output() ) then write(string1,*)trim(varname)//': replacing missing_value ',spvalR8,' with ',MISSING_R8 call error_handler(E_MSG,'get_var_3d',string1) endif io1 = nf90_get_att(ncid, VarID, 'scale_factor', scale_factor) io2 = nf90_get_att(ncid, VarID, 'add_offset' , add_offset) if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then where (var3d /= MISSING_R8) var3d = var3d * scale_factor + add_offset elseif (io1 == NF90_NOERR) then where (var3d /= MISSING_R8) var3d = var3d * scale_factor elseif (io2 == NF90_NOERR) then where (var3d /= MISSING_R8) var3d = var3d + add_offset endif else write(string1,*) trim(varname)//' has unsupported (by DART) xtype of', xtype call error_handler(E_ERR,'get_var_3d',string1,source,revision,revdate) endif end subroutine get_var_3d function get_model_time() type(time_type) :: get_model_time if ( .not. module_initialized ) call static_init_model get_model_time = model_time end function get_model_time function findVarIndex(varstring, caller) character(len=*), intent(in) :: varstring character(len=*), intent(in) :: caller integer :: findVarIndex integer :: i findVarIndex = -1 ! Skip to the right variable VARTYPES : do i = 1,nfields findVarIndex = i if ( trim(progvar(i)%varname) == varstring) exit VARTYPES enddo VARTYPES if (findVarIndex < 1) then write(string1,*) trim(caller)//' cannot find "'//trim(varstring)//'" in list of DART state variables.' call error_handler(E_ERR,'findVarIndex',string1,source,revision,revdate) endif end function findVarIndex subroutine update_water_table_depth( ivar, state_vector, ncFileID, filename, dart_time) !------------------------------------------------------------------ ! 'ivar' points to ZWT portion of the DART state vector, which has ! been updated by the assimilation at this point. ! ! ZWT is calculated from WA and WT ... so we have to update those ! CLM variables based on the new ZWT from the assimilation. ! Simply updating ZWT will have no effect because upon restart ! CLM will calculate ZWT given the same old WA and WT. ! ! So, this routine updates WA and WT based on the new ZWT. ! ! As defined in the CLM restart files: ! double WA(column) ; ! WA:long_name = "water in the unconfined aquifer" ; ! WA:units = "mm" ; ! double WT(column) ; ! WT:long_name = "total water storage" ; ! WT:units = "mm" ; ! double ZWT(column) ; ! ZWT:long_name = "water table depth" ; ! ZWT:units = "m" ; integer, intent(in) :: ivar real(r8), intent(in) :: state_vector(:) integer, intent(in) :: ncFileID character(len=*), intent(in) :: filename type(time_type), intent(in) :: dart_time ! temp space to hold data while we are writing it integer :: i, ni real(r8), allocatable, dimension(:) :: wa,wt,zwt real(r8), allocatable, dimension(:) :: zwt_pr !kyh11142014 integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs character(len=NF90_MAX_NAME) :: varname integer :: VarID, ncNdims, dimlen type(time_type) :: file_time ! FIXME ... check that dart_time and file_time agree ! TJH FIXME ... this algorithm is insufficient. Sean Swenson advised ! TJH FIXME ... that I should look into the code of the soil hydrology module ! TJH FIXME ... cesm1_1_1/models/lnd/clm/src/biogeophys/SoilHydrologyMod.F90 ! TJH FIXME ... for implementation. It gets more complicated when the water table ! TJH FIXME ... is more than the aquifer can hold and it starts wetting the soil. ! The simple test is to read in one file and write it out and compare ! to the original. No assimilation involved. Result: ! ! WA has min/max differences of -8.7606 0.00361754 (mm) ! WT has min/max differences of -131.581 1020.34 (mm) ! ! TJH FIXME ... So ... no dice. ! TJH FIXME ... Routine will error out until a better algorithm is in place. ! TJH FIXME ... AFTER TESTING, REMOVE THIS CHECK. write(string1,*)'WARNING: This routine is fundamentally untested.' write(string2,*)'WARNING: Exercise great care and check results thoroughly.' call error_handler(E_ERR, 'update_water_table_depth', string1, & source, revision, revdate, text2=string2) varname = trim(progvar(ivar)%varname) string2 = trim(filename)//' '//trim(varname) ! Ensure netCDF variable is conformable with progvar quantity. ! The TIME and Copy dimensions are intentionally not queried ! by looping over the dimensions stored in the progvar type. call nc_check(nf90_inq_varid(ncFileID, varname, VarID), & 'update_water_table_depth', 'inq_varid '//trim(string2)) call nc_check(nf90_inquire_variable(ncFileID,VarID,dimids=dimIDs,ndims=ncNdims), & 'update_water_table_depth', 'inquire '//trim(string2)) DimCheck : do i = 1,progvar(ivar)%numdims write(string1,'(''inquire dimension'',i2,A)') i,trim(string2) call nc_check(nf90_inquire_dimension(ncFileID, dimIDs(i), len=dimlen), & 'update_water_table_depth', string1) if ( dimlen /= progvar(ivar)%dimlens(i) ) then write(string1,*)trim(string2),' dim/dimlen ',i,dimlen,& ' not ',progvar(ivar)%dimlens(i) write(string2,*)' but it should be.' call error_handler(E_ERR, 'update_water_table_depth', string1, & source, revision, revdate, text2=string2) endif enddo DimCheck ni = progvar(ivar)%dimlens(1) ! number of CLM columns, in this context allocate(zwt(ni), wa(ni), wt(ni)) allocate(zwt_pr(ni)) !kyh11142014 ! Default values from ! cesm1_1_1/models/lnd/clm/src/main/mkarbinitMod.F90 ! FIXME ... should we read WA,WT originally and then update ... ! WA,WT have non-missing values in lake columns (the default ... 5000.0). wa(:) = 5000.0_r8 wt(:) = 5000.0_r8 zwt(:) = 0.0_r8 ! When called with a 4th argument, vector_to_prog_var() replaces the DART ! missing code with the value in the corresponding variable in the netCDF file. ! Any clamping to physically meaningful values occurrs in vector_to_prog_var. !call vector_to_prog_var(state_vector, ivar, zwt, ncFileID) !original call vector_to_prog_var(state_vector, ivar, zwt, zwt_pr, ncFileID) !kyh11142014 PARTITION: do i = 1,ni if (cols1d_ityplun(i) == LAKE) cycle PARTITION if ( zwt(i) == progvar(ivar)%spvalr8 ) then wa(i) = progvar(ivar)%spvalr8 wt(i) = progvar(ivar)%spvalr8 write(*,*)'FIXME hammering the value of wa, wt for column ',i elseif ( zwt(i) < 3.8019_r8 ) then wa(i) = 5000.0_r8 ! mm wt(i) = 5000.0_r8 + (3.8019_r8 - zwt(i))*400.0_r8 else wa(i) = (25.0_r8 - (zwt(i) - 3.8019_r8))*200.0_r8 wt(i) = wa(i) endif enddo PARTITION ! Stuff the updated ZWT values into the CLM netCDF file. call nc_check(nf90_put_var(ncFileID, VarID, zwt), & 'update_water_table_depth', 'put_var '//trim(varname)) ! Make note that the variable has been updated by DART ! call nc_check(nf90_Redef(ncFileID), & ! 'update_water_table_depth', 'redef '//trim(filename)) ! call nc_check(nf90_put_att(ncFileID, VarID,'DART','variable modified by DART'), & ! 'update_water_table_depth', 'modified '//trim(varname)) ! call nc_check(nf90_enddef(ncfileID), & ! 'update_water_table_depth','state enddef '//trim(filename)) ! Stuff the updated WA values into the CLM netCDF file. call nc_check(nf90_inq_varid(ncFileID, 'WA', VarID), & 'update_water_table_depth', 'inq_varid WA '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, wa), & 'update_water_table_depth', 'put_var WA '//trim(filename)) ! Make note that the variable has been updated by DART ! call nc_check(nf90_Redef(ncFileID), & ! 'update_water_table_depth', 'redef '//trim(filename)) ! call nc_check(nf90_put_att(ncFileID, VarID,'DART','variable modified by DART'), & ! 'update_water_table_depth', 'modified WA '//trim(filename)) ! call nc_check(nf90_enddef(ncfileID), & ! 'update_water_table_depth','state enddef '//trim(filename)) ! Stuff the updated WT values into the CLM netCDF file. call nc_check(nf90_inq_varid(ncFileID, 'WT', VarID), & 'update_water_table_depth', 'inq_varid WT '//trim(filename)) call nc_check(nf90_put_var(ncFileID, VarID, wt), & 'update_water_table_depth', 'put_var WT '//trim(filename)) ! Make note that the variable has been updated by DART ! call nc_check(nf90_Redef(ncFileID), & ! 'update_water_table_depth', 'redef '//trim(filename)) ! call nc_check(nf90_put_att(ncFileID, VarID,'DART','variable modified by DART'),& ! 'update_water_table_depth', 'modified WT '//trim(filename)) ! call nc_check(nf90_enddef(ncfileID), & ! 'update_water_table_depth','state enddef '//trim(filename)) deallocate(zwt, wa, wt) end subroutine update_water_table_depth subroutine get_brightness_temperature(state_time, location, metadata, obs_val, istatus) ! This is THE forward observation operator. Given the state and a location, return the value ! The parts of the state required for this forward operator are not required to ! be part of the DART state vector. They are currently directly harvested from the CLM ! restart file. As such, the posteriors are not informative. use radiative_transfer_mod, only : ss_model use rtm_main, only : rtm_main1 !kyh04022014 type(time_type), intent(in) :: state_time ! valid time of DART state type(location_type), intent(in) :: location ! observation location real(r8), dimension(:), intent(in) :: metadata real(r8), intent(out) :: obs_val ! model estimate of observation value integer, intent(out) :: istatus ! status of the calculation integer, parameter :: N_FREQ = 1 ! observations come in one frequency at a time integer, parameter :: N_POL = 2 ! code automatically computes both polarizations real(r8), parameter :: density_h2o = 1000.0_r8 ! Water density Kg/m3 ! variables required by ss_snow() routine real(r4), allocatable, dimension(:,:) :: y ! 2D array of snow properties !real(r4) :: aux_ins(5) ! [nsnowlyrs, ground_T, soilsat, poros, proportionality] !original real(r4) :: aux_ins(15) ! [nsnowlyrs, ground_T, soilsat, poros, proportionality,soil_liq,soil_ice, & ! sandf, clayf, Tcanopy, lai_vege, forc_pbot, forc_t, forc_rh, vf] !kyh04042014 integer :: ctrl(4) ! [n_lyrs, n_aux_ins, n_snow_ins, n_freq] real(r4) :: freq( N_FREQ) ! frequencies at which calculations are to be done real(r4) :: tetad(N_FREQ) ! incidence angle of satellite real(r4) :: tb_ubc(N_POL,N_FREQ) ! upper boundary condition brightness temperature real(r4) :: tb_out(N_POL,N_FREQ) ! calculated brightness temperature - output real(r4) :: stickiness_snow, b_prime_vege, x_lambda_vege !kyh11202014 real(r4) :: scatt_albedo_vege !kyh03052016 real(r8) :: tran_can_v_vege, tran_can_h_vege !kyh03092016 integer :: year_kyh, month_kyh, day_kyh, hour_kyh, minute_kyh, second_kyh !kyh04042014 ! support variables integer :: ncid character(len=256) :: filename integer, allocatable, dimension(:) :: columns_to_get real(r4), allocatable, dimension(:) :: tb real(r8), allocatable, dimension(:) :: weights real(r8), dimension(LocationDims) :: loc real(r8) :: loc_lon, loc_lat real(r8), allocatable, dimension(:) :: tran_can_v !kyh03092016 real(r8), allocatable, dimension(:) :: tran_can_h !kyh03092016 integer :: ilonmin(1), ilatmin(1) ! need to be array-valued for minloc intrinsic integer :: ilon, ilat, icol, ncols ! AMSR-E Tb observation metadata integer :: ens_index ! Ensemble member number integer :: landcovercode ! for this location - future use real(r8) :: frequency ! observation frequency real(r8) :: footprint ! observation footprint character :: polarization ! observation polarization real(r8) :: incidence_angle ! satellite incidence angle istatus = 1 obs_val = MISSING_R8 ! FIXME ... determine lon/lat indices ! Poor method ... will not work for single column case nor ! for irregular/unstructured grids might not work well over poles ... loc = get_location(location) loc_lon = loc(1) ! longitude of observation (in degrees) loc_lat = loc(2) ! latitude of observation (in degrees) ilonmin = minloc( abs(loc_lon-LON) ) ilatmin = minloc( abs(loc_lat-LAT) ) ilon = ilonmin(1) ilat = ilatmin(1) ncols = get_ncols_in_gridcell(ilon,ilat) ! Early return if there are no CLM columns at this location. ! Forward operator returns 'failed', but life goes on. if (ncols == 0) then if ((debug > 3) .and. do_output()) then write(string1, *) 'gridcell ilon/ilat (',ilon,ilat,') has no CLM columns - skipping.' write(string2, '(''obs lon,lat ('',f12.6,'','',f12.6,'')'')') loc_lon, loc_lat call error_handler(E_MSG,'get_brightness_temperature',string1,text2=string2) endif return endif !allocate( columns_to_get(ncols), tb(ncols), weights(ncols) ) !original allocate( columns_to_get(ncols), tb(ncols), weights(ncols), tran_can_v(ncols), tran_can_h(ncols) ) !kyh03092016 columns_to_get(:) = -1 tb(:) = 0.0_r4 weights(:) = 0.0_r8 tran_can_v(:) = 0.0_r8 !kyh03092016 tran_can_h(:) = 0.0_r8 !kyh03092016 call get_colids_in_gridcell(ilon, ilat, columns_to_get) ! FIXME Presently skipping gridcells with lakes. ! get_column_snow() must also modified to use bulk snow formulation for lakes. if ( any(cols1d_ityplun(columns_to_get) == LAKE)) then if ((debug > 1) .and. do_output()) then write(string1, *) 'gridcell ilon/ilat (',ilon,ilat,') has a lake - skipping.' write(string2, '(''obs lon,lat ('',f12.6,'','',f12.6,'')'')') loc_lon, loc_lat call error_handler(E_MSG,'get_brightness_temperature',string1,text2=string2) endif !deallocate(columns_to_get, tb, weights) !original deallocate(columns_to_get, tb, weights, tran_can_v, tran_can_h) !kyh03092016 return endif ! TJH FIXME - no bulletproofing on order ... landcovercode = int(metadata(1)) frequency = metadata(2) footprint = metadata(3) if ( metadata(4) > 0.0_r8 ) then polarization = 'H' else polarization = 'V' endif incidence_angle = metadata(5) ens_index = int(metadata(6)) if ((debug > 99) .and. do_output()) then write(*,*)'TJH debug ... computing gridcell ',ilon,ilat write(*,*)'TJH debug ... gridcell has columns ',columns_to_get write(*,*)'TJH debug ... ens_index ',ens_index write(*,*)'TJH debug ... landcovercode ',landcovercode write(*,*)'TJH debug ... frequency ',frequency write(*,*)'TJH debug ... footprint ',footprint write(*,*)'TJH debug ... polarization ',polarization write(*,*)'TJH debug ... incidence_angle ',incidence_angle endif tetad(:) = incidence_angle freq(:) = frequency ! need to know which restart file to use to harvest information call build_clm_instance_filename(ens_index, state_time, filename) call nc_check(nf90_open(trim(filename), NF90_NOWRITE, ncid), & 'get_brightness_temperature','open '//trim(filename)) ! Loop over all columns in the gridcell that has the right location. SNOWCOLS : do icol = 1,ncols weights(icol) = cols1d_wtxy( columns_to_get(icol)) ! relative weight of column call get_column_snow(ncid, filename, columns_to_get(icol)) ! allocates snowcolumn if ( (debug > 2) .and. do_output() ) then if (snowcolumn%nlayers < 1) then write(string1, *) 'column (',columns_to_get(icol),') has no snow' call error_handler(E_MSG,'get_brightness_temperature',string1) else write(*,*)'nprops ',snowcolumn%nprops write(*,*)'nlayers ',snowcolumn%nlayers write(*,*)'t_grnd ',snowcolumn%t_grnd write(*,*)'soilsat ',snowcolumn%soilsat write(*,*)'soilpor ',snowcolumn%soilporos write(*,*)'proconst ',snowcolumn%propconst write(*,*)'thickness',snowcolumn%thickness write(*,*)'density ',snowcolumn%density write(*,*)'radius ',snowcolumn%grain_radius write(*,*)'liqwater ',snowcolumn%liquid_water write(*,*)'temp ',snowcolumn%temperature write(*,*)'soil_liq ',snowcolumn%soil_liq !kyh04032014 write(*,*)'soil_ice ',snowcolumn%soil_ice !kyh04032014 write(*,*)'sandf ',snowcolumn%sandf !kyh04032014 write(*,*)'clayf ',snowcolumn%clayf !kyh04032014 write(*,*)'Tcanopy ',snowcolumn%Tcanopy !kyh04032014 write(*,*)'lai_vege ',snowcolumn%lai_vege !kyh04032014 write(*,*)'forc_pbot ',snowcolumn%forc_pbot !kyh04032014 write(*,*)'forc_t ',snowcolumn%forc_t !kyh04032014 write(*,*)'forc_rh ',snowcolumn%forc_rh !kyh04032014 write(*,*)'vf ',snowcolumn%vf !kyh04032014 write(*,*)'stickiness',snowcolumn%stickiness !kyh11202014 write(*,*)'b_prime',snowcolumn%b_prime !kyh11202014 write(*,*)'x_lambda',snowcolumn%x_lambda !kyh11202014 endif endif if ( snowcolumn%nlayers == 0 ) then ! If there is no snow, the ss_model will calculate the brightness ! temperature of the bare soil. To indicate this, aux_ins(1) must ! be 0 and ctrl(1) must be 1 ctrl(1) = 1 else ctrl(1) = snowcolumn%nlayers endif ctrl(2) = 0 ! not used as far as I can tell ctrl(3) = snowcolumn%nprops ctrl(4) = N_FREQ aux_ins(1) = real(snowcolumn%nlayers,r4) aux_ins(2) = snowcolumn%t_grnd aux_ins(3) = snowcolumn%soilsat aux_ins(4) = snowcolumn%soilporos aux_ins(5) = 0.5_r4 ! FIXME - hardwired !--------------------------------------------------------------kyh04032014 aux_ins(6) = snowcolumn%soil_liq aux_ins(7) = snowcolumn%soil_ice aux_ins(8) = snowcolumn%sandf aux_ins(9) = snowcolumn%clayf aux_ins(10) = snowcolumn%Tcanopy aux_ins(11) = snowcolumn%lai_vege aux_ins(12) = snowcolumn%forc_pbot aux_ins(13) = snowcolumn%forc_t aux_ins(14) = snowcolumn%forc_rh aux_ins(15) = snowcolumn%vf !--------------------------------------------------------------kyh04032014 !--------------------------------------------------------------RTM parameters !kyh11202014 stickiness_snow = snowcolumn%stickiness b_prime_vege = snowcolumn%b_prime x_lambda_vege = snowcolumn%x_lambda !--------------------------------------------------------------RTM parameters !kyh11202014 scatt_albedo_vege = snowcolumn%scatt_albedo !kyh03052016 !------------------------------------------------------------------- ! Ally's description of the y(:,4) variable. ! LWC in kg/m2 --> kg is mass of liquid water ! --> m2 is surface area ! LWC 'kg/m2' by water density 'kg/m3', we get depth ! of liquid water (m). Then, (depth of liquid water (m) / depth of snowpack (m)) ! provides the fraction of LWC (m water/m snowpack). ! Since, both liquid water and snowpack has same surface area (m2), ! we can also express it as 'm3/m3'. allocate( y(ctrl(1), snowcolumn%nprops) ) if ( aux_ins(1) > 0 ) then y(:,1) = snowcolumn%thickness y(:,2) = snowcolumn%density y(:,3) = snowcolumn%grain_radius * 2.0_r4 / 1000000.0_r4 ! need meters (from microns) y(:,4) = snowcolumn%liquid_water / (density_h2o * snowcolumn%thickness) y(:,5) = snowcolumn%temperature else ! dummy values for bare ground y(:,:) = 0.0_r4 endif ! FIXME Ally ... if you have a better way to specify/determine, ! tb_ubc do it here. Call your atmospheric model to calculate it. tb_ubc(:,N_FREQ) = (/ 2.7_r4, 2.7_r4 /) ! If the landcovercode indicates that you want to use ! a different radiative transfer model ... implement it here. ! this will involve changing the following 'if' statement. if (landcovercode >= 0 ) then ! the tb_out array contains the calculated brightness temperature outputs ! at each polarization (rows) and frequency (columns). !call ss_model(ctrl, freq, tetad, y, tb_ubc, aux_ins, tb_out) !original !-----------------------------------------------------------------------------kyh04042014 call get_date(state_time, year_kyh, month_kyh, day_kyh, hour_kyh, minute_kyh, second_kyh) !call rtm_main1(ctrl, freq, tetad, y, aux_ins, tb_out, month_kyh) !call rtm_main1(ctrl, freq, tetad, y, aux_ins, tb_out, month_kyh, & ! stickiness_snow, b_prime_vege, x_lambda_vege) !kyh11202014 !call rtm_main1(ctrl, freq, tetad, y, aux_ins, tb_out, month_kyh, & ! stickiness_snow, b_prime_vege, x_lambda_vege, scatt_albedo_vege) !kyh03052016 call rtm_main1(ctrl, freq, tetad, y, aux_ins, tb_out, month_kyh, & stickiness_snow, b_prime_vege, x_lambda_vege, scatt_albedo_vege, & tran_can_v_vege, tran_can_h_vege) !kyh03092016 !-----------------------------------------------------------------------------kyh04042014 else ! call to alternative radiative transfer model goes here. endif if ((debug > 2) .and. do_output()) then write(*,*)'column ', columns_to_get(icol),' tb_out is ',tb_out endif if (polarization == 'H') then tb(icol) = tb_out(1,1) ! second dimension is only 1 frequency tran_can_h(icol) = tran_can_h_vege !kyh03092016 else tb(icol) = tb_out(2,1) ! second dimension is only 1 frequency tran_can_v(icol) = tran_can_v_vege !kyh03092016 endif deallocate( y ) call destroy_column_snow() enddo SNOWCOLS call nc_check(nf90_close(ncid), 'get_brightness_temperature','close '//trim(filename)) ! FIXME ... account for heterogeneity somehow ... ! must aggregate all columns in the gridcell ! area-weight the average obs_val = sum(tb * weights) / sum(weights) tran_can_v_grd(ilon, ilat, ens_index) = sum(tran_can_v * weights) / sum(weights) !kyh03092016 tran_can_h_grd(ilon, ilat, ens_index) = sum(tran_can_h * weights) / sum(weights) !kyh03092016 !--------------------------------------------------------------------kyh_debug !write(*,*) 'ens_index=', ens_index !write(*,*) 'tran_can_v_grd=', tran_can_v_grd(ilon, ilat, ens_index) !--------------------------------------------------------------------kyh_debug if ((debug > 1) .and. do_output()) then write(*,*)'tb for all columns is ',tb write(*,*)'weights for all columns is ',weights write(*,*)'(weighted) obs value is ',obs_val write(*,*) endif !deallocate( columns_to_get, tb, weights ) !original deallocate( columns_to_get, tb, weights, tran_can_v, tran_can_h ) !kyh03092016 istatus = 0 end subroutine get_brightness_temperature subroutine get_column_snow(ncid, filename, snow_column ) ! Read all the variables needed for the radiative transfer model as applied ! to a single CLM column. ! ! The treatment of snow-related variables is complicated. ! The SNLSNO variable defines the number of snow layers with valid values. ! HOWEVER, if the snow depth is < 0.01 m, the snow is not represented by a layer, ! so the SNLSNO(i) is zero even though there is a trace of snow. ! Even a trace amount of snow results in some sort of snow cover fraction. ! ! Lakes are treated differently. ! The SNLSNO(i) is always zero, even though there is snow. ! The snow over lakes is wholly contained in the bulk formulation variables ! as opposed to the snow layer variables. ! float WATSAT(levgrnd, lat, lon) ; ! WATSAT:long_name = "saturated soil water content (porosity)" ; ! WATSAT:units = "mm3/mm3" ; ! WATSAT:_FillValue = 1.e+36f ; ! WATSAT:missing_value = 1.e+36f ; integer, intent(in) :: ncid character(len=*), intent(in) :: filename integer, intent(in) :: snow_column real(r8) :: t_grnd(1) ! ground temperature !------------------------------------------------------------kyh04042014 real(r4) :: sandfrac_c(1) ! sand fraction (fraction) real(r4) :: clayfrac_c(1) ! clay fraction (fraction) real(r4) :: t_veg_c(1) ! vegetation physical temperature (K) real(r4) :: elai_c(1) ! leaf area index real(r4) :: forc_pbot_c(1) ! ground level pressure (Pa) real(r4) :: forc_t_c(1) ! ground level air temperature (K) real(r4) :: forc_rh_c(1) ! relative humidity real(r4) :: vf_c(1) ! vegetated area fraction !------------------------------------------------------------kyh04042014 !------------------------------------------------------------RTM parameters !kyh11202014 real(r4) :: stickiness_c(1) !snowpack stickiness (-) real(r4) :: b_prime_c(1) !vegetation RTM coefficient (-) real(r4) :: x_lambda_c(1) !vegetation RTM coefficient (-) !------------------------------------------------------------RTM parameters !kyh11202014 real(r4) :: scatt_albedo_c(1) !vegetation single scattering albedo !kyh03052016 integer :: snlsno(1) ! number of snow layers real(r8), allocatable, dimension(:) :: h2osoi_liq, h2osoi_ice, t_soisno real(r8), allocatable, dimension(:) :: dzsno, zsno, zisno, snw_rds integer :: varid, ilayer, nlayers, ij integer, dimension(2) :: ncstart, nccount ! Get the (scalar) number of active snow layers for this column. call nc_check(nf90_inq_varid(ncid,'SNLSNO', varid), & 'get_column_snow', 'inq_varid SNLSNO'//trim(filename)) call nc_check(nf90_get_var( ncid, varid, snlsno, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var SNLSNO '//trim(filename)) nlayers = abs(snlsno(1)) ! Set some return values ! FIXME ... soilsat is a hardwired value ! FIXME ... soilporos is a hardwired value ! FIXME ... propconst is a hardwired value snowcolumn%nprops = 5 snowcolumn%nlayers = nlayers snowcolumn%soilsat = 0.3 ! aux_ins(3) soil saturation [fraction] snowcolumn%soilporos = 0.4 ! aux_ins(4) soil porosity [fraction] snowcolumn%propconst = 0.5 ! aux_ins(5) proportionality between grain size & correlation length. ! Get the ground temperature for this column. ! double T_GRND(column); long_name = "ground temperature" ; units = "K" ; call nc_check(nf90_inq_varid(ncid,'T_GRND', varid), & 'get_column_snow', 'inq_varid T_GRND '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, t_grnd, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var T_GRND '//trim(filename)) !--------------------------------------------------kyh04042014 ! Get other information for this column ! double SANDFRAC_C(column); long_name='sand fraction in a column', units='fraction') call nc_check(nf90_inq_varid(ncid,'SANDFRAC_C', varid), & 'get_column_snow', 'inq_varid SANDFRAC_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, sandfrac_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var SANDFRAC_C '//trim(filename)) ! double CLAYFRAC_C(column); long_name='clay fraction in a column', units='fraction') call nc_check(nf90_inq_varid(ncid,'CLAYFRAC_C', varid), & 'get_column_snow', 'inq_varid CLAYFRAC_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, clayfrac_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var CLAYFRAC_C '//trim(filename)) ! double T_VEG_C(column); long_name='vegetation physical temperature in a column', units='K') call nc_check(nf90_inq_varid(ncid,'T_VEG_C', varid), & 'get_column_snow', 'inq_varid T_VEG_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, t_veg_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var T_VEG_C '//trim(filename)) ! double ELAI_C(column); long_name='leaf area index in a column', units='-') call nc_check(nf90_inq_varid(ncid,'ELAI_C', varid), & 'get_column_snow', 'inq_varid ELAI_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, elai_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var ELAI_C '//trim(filename)) ! double VF_C(column); long_name='vegetated area fraction in a column', units='fraction') call nc_check(nf90_inq_varid(ncid,'VF_C', varid), & 'get_column_snow', 'inq_varid VF_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, vf_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var VF_C '//trim(filename)) ! double FORC_PBOT_C(column); long_name='atmospheric pressure', units='Pa') call nc_check(nf90_inq_varid(ncid,'FORC_PBOT_C', varid), & 'get_column_snow', 'inq_varid FORC_PBOT_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, forc_pbot_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var FORC_PBOT_C '//trim(filename)) ! double FORC_T_C(column); long_name='atmospheric temperature', units='K') call nc_check(nf90_inq_varid(ncid,'FORC_T_C', varid), & 'get_column_snow', 'inq_varid FORC_T_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, forc_t_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var FORC_T_C '//trim(filename)) ! double FORC_RH_C(column); long_name='atmospheric relative humidity', units='%') call nc_check(nf90_inq_varid(ncid,'FORC_RH_C', varid), & 'get_column_snow', 'inq_varid FORC_RH_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, forc_rh_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var FORC_RH_C '//trim(filename)) !--------------------------------------------------kyh04042014 !--------------------------------------------------------------------------------------------------kyh12182014 (RTM parameters) ! double sTICKINESS_C(column); long_name=''Snowpack stickiness', units='unitless') call nc_check(nf90_inq_varid(ncid,'STICKINESS_C', varid), & 'get_column_snow', 'inq_varid STICKINESS_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, stickiness_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var STICKINESS_C '//trim(filename)) ! double B_PRIME_C(column); long_name='Vegetation RTM coefficient', units='unitless') call nc_check(nf90_inq_varid(ncid,'B_PRIME_C', varid), & 'get_column_snow', 'inq_varid B_PRIME_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, b_prime_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var B_PRIME_C '//trim(filename)) ! double X_LAMBDA_C(column); long_name='Vegetation RTM coefficient', units='unitless') call nc_check(nf90_inq_varid(ncid,'X_LAMBDA_C', varid), & 'get_column_snow', 'inq_varid X_LAMBDA_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, x_lambda_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var X_LAMBDA_C '//trim(filename)) !--------------------------------------------------------------------------------------------------kyh12182014 (RTM parameters) !-------------------------------------------------------------------------------------------------kyh03052016 ! double SCATT_ALBEDO_C(column); long_name='Vegetation single scattering albedo', units='unitless') call nc_check(nf90_inq_varid(ncid,'SCATT_ALBEDO_C', varid), & 'get_column_snow', 'inq_varid SCATT_ALBEDO_C '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, scatt_albedo_c, start=(/ snow_column /), count=(/ 1 /)), & 'get_column_snow', 'get_var SCATT_ALBEDO_C '//trim(filename)) !-------------------------------------------------------------------------------------------------kyh03052016 ! FIXME ... lake columns use a bulk formula for snow if (cols1d_ityplun(snow_column) == LAKE ) return ! we are a lake ! double H2OSOI_LIQ(column, levtot); long_name = "liquid water" ; units = "kg/m2" ; ! double H2OSOI_ICE(column, levtot); long_name = "ice lens" ; units = "kg/m2" ; ! double T_SOISNO( column, levtot); long_name = "soil-snow temperature" ; units = "K" ; allocate(h2osoi_liq(nlevtot), h2osoi_ice(nlevtot), t_soisno(nlevtot)) ncstart = (/ 1, snow_column /) nccount = (/ nlevtot, 1 /) call nc_check(nf90_inq_varid(ncid,'T_SOISNO', varid), & 'get_column_snow', 'inq_varid T_SOISNO '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, t_soisno, start=ncstart, count=nccount), & 'get_column_snow', 'get_var T_SOISNO '//trim(filename)) call nc_check(nf90_inq_varid(ncid,'H2OSOI_LIQ', varid), & 'get_column_snow', 'inq_varid H2OSOI_LIQ '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, h2osoi_liq, start=ncstart, count=nccount), & 'get_column_snow', 'get_var H2OSOI_LIQ '//trim(filename)) call nc_check(nf90_inq_varid(ncid,'H2OSOI_ICE', varid), & 'get_column_snow', 'inq_varid H2OSOI_ICE '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, h2osoi_ice, start=ncstart, count=nccount), & 'get_column_snow', 'get_var H2OSOI_ICE '//trim(filename)) ! double DZSNO(column, levsno); long_name = "snow layer thickness" ; units = "m" ; ! double ZSNO(column, levsno); long_name = "snow layer depth" ; units = "m" ; ! double ZISNO(column, levsno); long_name = "snow interface depth" ; units = "m" ; ! double snw_rds(column, levsno); long_name = "snow layer effective radius" ; units = "um" ; allocate(dzsno(nlevsno), zsno(nlevsno), zisno(nlevsno), snw_rds(nlevsno)) ncstart = (/ 1, snow_column /) nccount = (/ nlevsno, 1 /) call nc_check(nf90_inq_varid(ncid,'DZSNO', varid), & 'get_column_snow', 'inq_varid DZSNO '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, dzsno, start=ncstart, count=nccount), & 'get_column_snow', 'get_var DZSNO '//trim(filename)) call nc_check(nf90_inq_varid(ncid,'ZSNO', varid), & 'get_column_snow', 'inq_varid ZSNO '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, zsno, start=ncstart, count=nccount), & 'get_column_snow', 'get_var ZSNO '//trim(filename)) call nc_check(nf90_inq_varid(ncid,'ZISNO', varid), & 'get_column_snow', 'inq_varid ZISNO '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, zisno, start=ncstart, count=nccount), & 'get_column_snow', 'get_var ZISNO '//trim(filename)) call nc_check(nf90_inq_varid(ncid,'snw_rds', varid), & 'get_column_snow', 'inq_varid snw_rds '//trim(filename)) call nc_check(nf90_get_var( ncid, varid, snw_rds, start=ncstart, count=nccount), & 'get_column_snow', 'get_var snw_rds '//trim(filename)) ! Print a summary so far if ((debug > 3) .and. do_output()) then write(*,*)'get_column_snow: raw CLM data for column ',snow_column write(*,*)' # of snow layers, column ityp, ground temp :', snlsno, cols1d_ityplun(snow_column), t_grnd write(*,*)' h2osoi_liq :', h2osoi_liq(1:nlevsno) write(*,*)' h2osoi_ice :', h2osoi_ice(1:nlevsno) write(*,*)' t_soisno :', t_soisno(1:nlevsno) write(*,*)' dzsno :', dzsno(1:nlevsno) write(*,*)' zsno :', zsno(1:nlevsno) write(*,*)' zisno :', zisno(1:nlevsno) write(*,*)' snw_rds :', snw_rds(1:nlevsno) endif allocate( snowcolumn%thickness(nlayers) , & snowcolumn%density(nlayers) , & snowcolumn%grain_radius(nlayers) , & snowcolumn%liquid_water(nlayers) , & snowcolumn%temperature(nlayers) ) ! Fill the output array ... finally snowcolumn%t_grnd = t_grnd(1) !-------------------------------------kyh04032014 snowcolumn%soil_liq = h2osoi_liq(6) !soil layer index starts from 6 (since snow layer index is from 1 to 5) snowcolumn%soil_ice = h2osoi_ice(6) snowcolumn%sandf = sandfrac_c(1) snowcolumn%clayf = clayfrac_c(1) snowcolumn%Tcanopy = t_veg_c(1) snowcolumn%lai_vege = elai_c(1) snowcolumn%forc_pbot= forc_pbot_c(1) snowcolumn%forc_t = forc_t_c(1) snowcolumn%forc_rh = forc_rh_c(1) snowcolumn%vf = vf_c(1) !-------------------------------------kyh04032014 !-------------------------------------RTM parameter !kyh11202014 snowcolumn%stickiness = stickiness_c(1) snowcolumn%b_prime = b_prime_c(1) snowcolumn%x_lambda = x_lambda_c(1) !-------------------------------------RTM parameter !kyh11202014 snowcolumn%scatt_albedo = scatt_albedo_c(1) !kyh03052016 ij = 0 do ilayer = (nlevsno-nlayers+1),nlevsno ij = ij + 1 snowcolumn%thickness(ij) = dzsno(ilayer) snowcolumn%density(ij) = (h2osoi_liq(ilayer) + h2osoi_ice(ilayer)) / dzsno(ilayer) snowcolumn%grain_radius(ij) = snw_rds(ilayer) snowcolumn%liquid_water(ij) = h2osoi_liq(ilayer) snowcolumn%temperature(ij) = t_soisno(ilayer) if ((debug > 3) .and. do_output()) & write(*,*)' get_column_snow: filling column ',snow_column, & ' layer ',ij,' with info from ilayer ',ilayer enddo deallocate(h2osoi_liq, h2osoi_ice, t_soisno) deallocate(dzsno, zsno, zisno, snw_rds) end subroutine get_column_snow subroutine destroy_column_snow if (associated(snowcolumn%thickness)) deallocate(snowcolumn%thickness) if (associated(snowcolumn%density)) deallocate(snowcolumn%density) if (associated(snowcolumn%grain_radius)) deallocate(snowcolumn%grain_radius) if (associated(snowcolumn%liquid_water)) deallocate(snowcolumn%liquid_water) if (associated(snowcolumn%temperature)) deallocate(snowcolumn%temperature) snowcolumn%nprops = 0 snowcolumn%nlayers = 0 snowcolumn%t_grnd = 0.0_r4 snowcolumn%soilsat = 0.0_r4 snowcolumn%soilporos = 0.0_r4 snowcolumn%propconst = 0.0_r4 !--------------------------------kyh04032014 snowcolumn%soil_liq = 0.0_r4 snowcolumn%soil_ice = 0.0_r4 snowcolumn%sandf = 0.0_r4 snowcolumn%clayf = 0.0_r4 snowcolumn%Tcanopy = 0.0_r4 snowcolumn%lai_vege = 0.0_r4 snowcolumn%forc_pbot = 0.0_r4 snowcolumn%forc_t = 0.0_r4 snowcolumn%forc_rh = 0.0_r4 snowcolumn%vf = 0.0_r4 !--------------------------------kyh04032014 end subroutine destroy_column_snow !====================================================================== subroutine build_clm_instance_filename(instance, state_time, filename) ! If the instance is 1, it could be a perfect model scenario ! or it could be the first instance of many. CESM has a different ! naming scheme for these. integer, intent(in) :: instance type(time_type), intent(in) :: state_time character(len=*), intent(out) :: filename integer :: year, month, day, hour, minute, second 100 format (A,'.clm2_',I4.4,'.r.',I4.4,'-',I2.2,'-',I2.2,'-',I5.5,'.nc') 110 format (A,'.clm2' ,'.r.',I4.4,'-',I2.2,'-',I2.2,'-',I5.5,'.nc') call get_date(state_time, year, month, day, hour, minute, second) second = second + minute*60 + hour*3600 write(filename,110) trim(casename),year,month,day,second if( file_exist(filename) ) then ! perfect model scenario if ( (debug > 99) .and. do_output()) then write(string1,*)'Running in a perfect model configuration with ',trim(filename) call error_handler(E_MSG, 'model_mod:build_clm_instance_filename', string1) endif else ! multi-instance scenario write(filename,100) trim(casename),instance,year,month,day,second if( .not. file_exist(filename) ) then write(string1,*)'Unable to create viable CLM restart filename:' call error_handler(E_ERR, 'model_mod:build_clm_instance_filename', & string1, text2=trim(filename)) endif endif end subroutine build_clm_instance_filename subroutine SetLocatorArrays() ! This function will create the relational table that will indicate how many ! and which columns pertain to the gridcells. A companion function will ! return the column indices that are needed to recreate the gridcell value. ! ! This fills the gridCellInfo(:,:) structure. ! given a gridcell, the gridCellInfo(:,:) structure will indicate how many and ! which columns are part of the gridcell. integer :: ilon, ilat, ij, iunit integer :: icol, currenticol(nlon,nlat) integer :: ipft, currentipft(nlon,nlat) gridCellInfo(:,:)%ncols = 0 gridCellInfo(:,:)%npfts = 0 ! Count up how many columns are in each gridcell do ij = 1,Ncolumn ilon = cols1d_ixy(ij) ilat = cols1d_jxy(ij) gridCellInfo(ilon,ilat)%ncols = gridCellInfo(ilon,ilat)%ncols + 1 enddo ! Count up how many pfts are in each gridcell do ij = 1,Npft ilon = pfts1d_ixy(ij) ilat = pfts1d_jxy(ij) gridCellInfo(ilon,ilat)%npfts = gridCellInfo(ilon,ilat)%npfts + 1 enddo ! Create storage for the list of column,pft indices do ilon = 1,nlon do ilat = 1,nlat if ( gridCellInfo(ilon,ilat)%ncols > 0 ) then allocate( gridCellInfo(ilon,ilat)%columnids( gridCellInfo(ilon,ilat)%ncols )) endif if ( gridCellInfo(ilon,ilat)%npfts > 0 ) then allocate( gridCellInfo(ilon,ilat)%pftids( gridCellInfo(ilon,ilat)%npfts )) endif enddo enddo ! Fill the column pointer arrays currenticol(:,:) = 0 do ij = 1,Ncolumn ilon = cols1d_ixy(ij) ilat = cols1d_jxy(ij) currenticol(ilon,ilat) = currenticol(ilon,ilat) + 1 icol = currenticol(ilon,ilat) if ( icol <= gridCellInfo(ilon,ilat)%ncols ) then gridCellInfo(ilon,ilat)%columnids(icol) = ij else write(string1,'(''gridcell('',i4,'','',i4,'') has at most '',i4,'' columns.'')') & ilon, ilat, gridCellInfo(ilon,ilat)%ncols write(string2,'(''Found '',i8,'' at dart index '',i12)') icol, ij call error_handler(E_ERR, 'SetLocatorArrays', string1, & source, revision, revdate, text2=string2) endif enddo ! Fill the pft pointer arrays currentipft(:,:) = 0 do ij = 1,Npft ilon = pfts1d_ixy(ij) ilat = pfts1d_jxy(ij) currentipft(ilon,ilat) = currentipft(ilon,ilat) + 1 ipft = currentipft(ilon,ilat) if ( ipft <= gridCellInfo(ilon,ilat)%npfts ) then gridCellInfo(ilon,ilat)%pftids(ipft) = ij else write(string1,'(''gridcell('',i4,'','',i4,'') has at most '',i4,'' pfts.'')') & ilon, ilat, gridCellInfo(ilon,ilat)%npfts write(string2,'(''Found '',i8,'' at dart index '',i12)') ipft, ij call error_handler(E_ERR, 'SetLocatorArrays', string1, & source, revision, revdate, text2=string2) endif enddo ! Check block if ((debug > 99) .and. do_output()) then iunit = open_file('gridcell_column_table.txt',form='formatted',action='write') do ilon = 1,nlon do ilat = 1,nlat if (gridCellInfo(ilon,ilat)%ncols > 0) then write(iunit,'(''gridcell'',i8,1x,i8,'' has '', i6, '' columns:'')') & ilon,ilat,gridCellInfo(ilon,ilat)%ncols write(iunit,*)gridCellInfo(ilon,ilat)%columnids endif enddo enddo call close_file(iunit) iunit = open_file('gridcell_pft_table.txt',form='formatted',action='write') do ilon = 1,nlon do ilat = 1,nlat if (gridCellInfo(ilon,ilat)%npfts > 0) then write(iunit,'(''gridcell'',i8,1x,i8,'' has '', i6, '' pfts : '')') & ilon,ilat,gridCellInfo(ilon,ilat)%npfts write(iunit,*)gridCellInfo(ilon,ilat)%pftids endif enddo enddo call close_file(iunit) endif end subroutine SetLocatorArrays !> SetVariableAttributes() converts the information in the variable_table !> to the progvar structure for each variable. !> If the numerical limit does not apply, it is set to MISSING_R8, even if !> it is the maximum that does not apply. subroutine SetVariableAttributes(ivar) integer, intent(in) :: ivar integer :: ios real(r8) :: minvalue, maxvalue progvar(ivar)%varname = trim(variable_table(ivar,VT_VARNAMEINDX)) progvar(ivar)%kind_string = trim(variable_table(ivar,VT_KINDINDX)) progvar(ivar)%dart_kind = get_raw_obs_kind_index( progvar(ivar)%kind_string ) progvar(ivar)%maxlevels = 0 progvar(ivar)%dimlens = 0 progvar(ivar)%dimnames = ' ' progvar(ivar)%spvalINT = -9999 ! from CESM/CLM clm_varcon.F90 progvar(ivar)%spvalR4 = 1.e36_r4 ! from CESM/CLM clm_varcon.F90 progvar(ivar)%spvalR8 = 1.e36_r8 ! from CESM/CLM clm_varcon.F90 progvar(ivar)%missingINT = MISSING_I progvar(ivar)%missingR4 = MISSING_R4 progvar(ivar)%missingR8 = MISSING_R8 progvar(ivar)%rangeRestricted = BOUNDED_NONE progvar(ivar)%minvalue = MISSING_R8 progvar(ivar)%maxvalue = MISSING_R8 progvar(ivar)%has_fill_value = .true. progvar(ivar)%has_missing_value = .true. progvar(ivar)%update = .false. if (variable_table(ivar,VT_ORIGININDX) == 'VECTOR') then progvar(ivar)%origin = trim(clm_vector_history_filename) elseif (variable_table(ivar,VT_ORIGININDX) == 'LAIIN') then progvar(ivar)%origin = trim(lai_nc) !kyh01192016 elseif (variable_table(ivar,VT_ORIGININDX) == 'HISTORY') then progvar(ivar)%origin = trim(clm_history_filename) else variable_table(ivar,VT_ORIGININDX) = 'RESTART' progvar(ivar)%origin = trim(clm_restart_filename) endif if ((variable_table(ivar,VT_STATEINDX) == 'UPDATE') .and. & (variable_table(ivar,VT_ORIGININDX) == 'RESTART')) progvar(ivar)%update = .true. ! set the default values minvalue = MISSING_R8 maxvalue = MISSING_R8 progvar(ivar)%minvalue = MISSING_R8 progvar(ivar)%maxvalue = MISSING_R8 ! If the character string can be interpreted as an r8, great. ! If not, there is no value to be used. read(variable_table(ivar,VT_MINVALINDX),*,iostat=ios) minvalue if (ios == 0) progvar(ivar)%minvalue = minvalue read(variable_table(ivar,VT_MAXVALINDX),*,iostat=ios) maxvalue if (ios == 0) progvar(ivar)%maxvalue = maxvalue ! rangeRestricted == BOUNDED_NONE == 0 ... unlimited range ! rangeRestricted == BOUNDED_BELOW == 1 ... minimum, but no maximum ! rangeRestricted == BOUNDED_ABOVE == 2 ... maximum, but no minimum ! rangeRestricted == BOUNDED_BOTH == 3 ... minimum and maximum if ( (progvar(ivar)%minvalue /= MISSING_R8) .and. & (progvar(ivar)%maxvalue /= MISSING_R8) ) then progvar(ivar)%rangeRestricted = BOUNDED_BOTH elseif (progvar(ivar)%maxvalue /= MISSING_R8) then progvar(ivar)%rangeRestricted = BOUNDED_ABOVE elseif (progvar(ivar)%minvalue /= MISSING_R8) then progvar(ivar)%rangeRestricted = BOUNDED_BELOW else progvar(ivar)%rangeRestricted = BOUNDED_NONE endif ! Check to make sure min is less than max if both are specified. if ( progvar(ivar)%rangeRestricted == BOUNDED_BOTH ) then if (maxvalue < minvalue) then write(string1,*)'&model_nml state_variable input error for ',trim(progvar(ivar)%varname) write(string2,*)'minimum value (',minvalue,') must be less than ' write(string3,*)'maximum value (',maxvalue,')' call error_handler(E_ERR,'SetVariableAttributes',string1, & source,revision,revdate,text2=trim(string2),text3=trim(string3)) endif endif end subroutine SetVariableAttributes !> FindDesiredTimeIndx() returns the index into the time array that matches !> the model_time from the CLM restart file. function FindDesiredTimeIndx(ncid, ntimes, varname) integer, intent(in) :: ncid integer, intent(in) :: ntimes character(len=*), intent(in) :: varname integer :: FindDesiredTimeIndx integer :: VarID real(r8), dimension(ntimes) :: mytimes character(len=NF90_MAX_NAME) :: attvalue type(time_type) :: thistime integer :: ios, itime, basedays, baseseconds integer :: iyear, imonth, iday, ihour, imin, isec FindDesiredTimeIndx = MISSING_I ! initialize to failure setting call nc_check(nf90_inq_varid(ncid, 'time', VarID), & 'FindDesiredTimeIndx:', 'inq_varid time '//varname) call nc_check(nf90_get_var( ncid, VarID, mytimes), & 'FindDesiredTimeIndx:', 'get_var time '//varname) call nc_check(nf90_get_att(ncid, VarID, 'units', attvalue), & 'FindDesiredTimeIndx:', 'time get_att units '//varname) ! time:units = "days since 2004-01-01 00:00:00" ; ! 1234567890 if (attvalue(1:10) /= 'days since') then write(string1,*)'expecting time units of [days since ... ]' write(string2,*)'read time units of ['//trim(attvalue)//']' call error_handler(E_ERR, 'FindDesiredTimeIndx:', string1, & source, revision, revdate, text2=string2) endif read(attvalue,'(11x,i4,5(1x,i2))',iostat=ios)iyear,imonth,iday,ihour,imin,isec if (ios /= 0) then write(string1,*)'Unable to read time units. Error status was ',ios write(string2,*)'expected "days since YYYY-MM-DD HH:MM:SS"' write(string3,*)'was "'//trim(attvalue)//'"' call error_handler(E_ERR, 'FindDesiredTimeIndx:', string1, & source, revision, revdate, text2=string2, text3=string3) endif thistime = set_date(iyear, imonth, iday, ihour, imin, isec) call get_time(thistime, baseseconds, basedays) ! convert each time to a DART time and compare to desired TIMELOOP : do itime = 1,ntimes iday = int(mytimes(itime)) isec = (mytimes(itime) - iday)*86400 thistime = set_time(baseseconds+isec, basedays+iday) if (thistime == model_time) then FindDesiredTimeIndx = itime exit TIMELOOP endif enddo TIMELOOP ! FIXME ... do we actually need a perfect match ... or do we just use the last ! one. History files may not have a perfect match to the restart file. if ( FindDesiredTimeIndx == MISSING_I ) then call print_time(model_time,str='model time is ',iunit=logfileunit) call print_time(model_time,str='model time is ') call print_date(model_time,str='model date is ',iunit=logfileunit) call print_date(model_time,str='model date is ') write(string1,*)'No time matching model_time found for '//trim(varname) call error_handler(E_ERR, 'FindDesiredTimeIndx:', string1, & source, revision, revdate ) endif if ((debug > 0) .and. do_output()) then write(string1,*)trim(varname)//' matching time index is ',FindDesiredTimeIndx call error_handler(E_MSG, 'FindDesiredTimeIndx:', string1) endif end function FindDesiredTimeIndx !=================================================================== ! End of model_mod !=================================================================== end module model_mod ! ! $URL: https://proxy.subversion.ucar.edu/DAReS/DART/branches/Tb/models/clm/model_mod.f90 $ ! $Id: model_mod.f90 7202 2014-10-08 15:29:19Z thoar $ ! $Revision: 7202 $ ! $Date: 2014-10-08 10:29:19 -0500 (Wed, 08 Oct 2014) $