diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 5607ee5e3d..2d581f9b00 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -6,4 +6,4 @@ Marco VAN HULTEN Ana GONZALEZ-NICOLAS Daniel CAVIEDES-VOULLIÈME Olga DOMBROWSKI - +Yorck EWERDWALBESLOH diff --git a/docs/_toc.yml b/docs/_toc.yml index 7f74453e61..659415e6bb 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -9,6 +9,8 @@ parts: title: Installing eCLM - file: users_guide/introduction/introduction title: Scientific Background + - file: users_guide/introduction/perturbation + title: Perturbation routines - caption: User's Guide chapters: diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md new file mode 100644 index 0000000000..9fd8ec73b5 --- /dev/null +++ b/docs/users_guide/introduction/perturbation.md @@ -0,0 +1,472 @@ +# Perturbation routines # + +This section describes **perturbation capabilities for eCLM-PDAF**. + +The implementation of the perturbation routines was first introduced +by Yorck Ewerdwalbesloh in +https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5. + +## Description + +eCLM implements perturbations for atmospheric forcing variables +(temperature, precipitation, longwave radiation, and shortwave +radiation) to support ensemble-based data assimilation in +eCLM-PDAF. The perturbation scheme maintains physical consistency +through cross-variable correlations and spatiotemporal coherence. + +### Physical Consistency + +Perturbed variables are cross-correlated to preserve physical +relationships between forcing fields. For example, a positive +perturbation of incoming shortwave radiation corresponds to a negative +perturbation of longwave radiation and a positive perturbation of +temperature. The cross-correlation structure follows Reichle et +al. (2007) and Han et al. (2014). + +### Spatial Correlation + +Perturbations are spatially correlated using an isotropic correlation +function based on grid cell separation distance. Precipitation uses a +shorter correlation length scale than other variables to reflect its +more localized spatial structure. Correlation parameters are +configurable. + +### Temporal Correlation + +Temporal correlation is generated using a first-order autoregressive +(AR-1) model following Evensen (2009). The default temporal +persistence corresponds to a decorrelation timescale of one day, which +is adjustable based on application requirements. + +### Implementation + +Perturbation fields are read from preprocessed noise files and applied +to atmospheric forcing variables within the eCLM model during runtime. + +## Key Components ## + +### 1. Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) ### + +Allows reading perturbed hydraulic parameters from input files instead +of computing them via pedotransfer functions. + +**Parameters:** +- `THETAS` - saturated water content +- `SHAPE_PARAM` - Brooks-Corey shape parameter (`bsw`) +- `PSIS_SAT` - saturated matric potential (`sucsat`) +- `KSAT` - saturated hydraulic conductivity (`xksat`) + +**Implementation:** +- Hydraulic parameters read from surface file +- Reads EITHER original parameters (for `nlevsoifl=10` soil layers) + and applies organic matter mixing OR: Reads adjusted parameters with + `_adj` suffix (for all `nlevgrnd` layers) and **overwrites** + parameters from organic matter mixing. +- Falls back to pedotransfer functions if parameters aren't in the file +- Modifies organic matter mixing to preserve perturbed parameter + values + +#### Surface File + +Soil hydraulic parameters are read from the **surface dataset file**, +which is specified in the `lnd_in` namelist file under key `fsurdat`. + +##### Sufrace File: Parameter Structure and Dimensions + +All soil hydraulic parameters are stored as three-dimensional arrays +with the structure: + +``` + (nlevgrnd, lsmlat, lsmlon) +``` + +Where: +- `nlevgrnd`: Number of soil layers (vertical dimension, typically 25 + layers in eCLM simulations, see namelist file key + `soil_layerstruct`) +- `lsmlat`: Number of latitude grid points +- `lsmlon`: Number of longitude grid points + +This 3D structure enables depth-varying, spatially-distributed soil +properties across the model domain. + +##### Surface File: Input Parameter Format + +Two options +- original parameters that undergo organic mixing +- adjusted parameters that are used by input values + +###### Original Parameters + +The surface file must contain the following baseline parameters: + +``` text + THETAS(nlevgrnd, lsmlat, lsmlon) + long_name: "Porosity" + units: "vol/vol" + + SHAPE_PARAM(nlevgrnd, lsmlat, lsmlon) + long_name: "Shape (b) parameter" + units: "unitless" + + PSIS_SAT(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated soil matric potential" + units: "mmH2O" + + KSAT(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated hydraulic conductivity" + units: "mm/s" +``` + +These parameters apply to the first 10 soil layers (`nlevsoifl=10`) and +undergo organic matter mixing. + +###### Adjusted Parameters for Data Assimilation + +When running with perturbed soil parameters (e.g., from ensemble data +assimilation), the surface file should contain adjusted parameters +with the `_adj` suffix: + +``` text + THETAS_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Porosity" + units: "vol/vol" + _FillValue: 1.e+30 + + SHAPE_PARAM_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Shape (b) parameter" + units: "unitless" + _FillValue: 1.e+30 + + PSIS_SAT_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated soil matric potential" + units: "mmH2O" + _FillValue: 1000.0 + + KSAT_adj(nlevgrnd, lsmlat, lsmlon) + long_name: "Saturated hydraulic conductivity" + units: "mm/s" + _FillValue: 1.e+30 +``` + +**Important Notes:** +- Adjusted parameters apply to **all** `nlevgrnd` layers +- When present, adjusted parameters **overwrite** the results of organic + matter mixing +- `PSIS_SAT_adj` uses a different fill value (1000.0) compared to other + parameters (1.e+30), reflecting special handling for undefined matric + potential values +- If adjusted parameters are not present in the surface file, eCLM falls + back to the original parameters or pedotransfer functions + +#### Parameter Priority and Fallback Option + +The code in `SoilStateInitTimeConstMod.F90` follows this priority order: + +1. **First priority**: Read `_adj` parameters (if present) → applies + to all `nlevgrnd` layers → **overwrites** organic matter mixing + results +2. **Second priority**: Read original parameters → applies to first 10 + layers (`nlevsoifl`) → undergoes organic matter mixing +3. **Fallback**: Use pedotransfer functions if parameters are absent + from the file + +This hierarchical approach ensures maximum flexibility for both +standard simulations and data assimilation applications. + +#### Note about the Brooks-Corey Shape Parameter #### + +When perturbed soil parameters are read from input files, the organic +matter mixing for `bsw` uses the file-read value instead of the +hard-coded Cosby et al. (1984) Table 5 formula (`2.91 + 0.159*clay`). + +### 2. Noise-Based Forcing Perturbation ### + +When having an ensemble run, the memory consumption is too large if +the forcing data is perturbed one by one so that I get a file for each +month for each member idea: perturb the forcings in the CLM sourcecode +with a noise file. for each perturbed variable (temperature, +precipitation, longwave and shortwave radiation), a stream is +introduced. + +Adds spatiotemporal noise to atmospheric forcing data for ensemble +data assimilation. + +**Modified files:** +- `shr_stream_mod.F90` +- `shr_strdata_mod.F90` +- `shr_dmodel_mod.F90` + +**Key Features:** +- Stores ensemble metadata in stream structure: + - `caseId` - ensemble member ID + - `dt` - forcing time resolution + - `numEns` - ensemble size that the perturbation file was created + for. Example: Running with 50 ensemble members with a noise file + created for an ensemble of up to 200 members. Then, `numEns` + should be set to 200 in the stream file for right usage of + perturbation dimension. +- **Time-shifting mechanism**: Different ensemble members read different temporal frames from the same noise file + - Formula: `frame = nt + caseId * 24/dt` + - Example: For 3-hourly data (`dt=3`), member 0 starts at frame + `nt`, member 1 at `nt+8`, member 2 at `nt+16`, etc. +- Detects noise fields by checking if model field name contains `"noise"` +- Adjusts time coordinate reading to account for ensemble-extended noise files + +#### Preparing Atmospheric Forcing Noise #### + +This section describes the workflow for generating and configuring +spatiotemporally correlated noise fields for atmospheric forcing +perturbations. + +##### Noise Generation Tool ##### + +The `correlatedNoise.py` +(https://github.com/HPSCTerrSys/eCLM_atm-forcing-generator/pull/8) +script generates monthly NetCDF files containing spatiotemporally +correlated noise fields for all ensemble members. + +**Key Configuration Parameters:** + +- **Ensemble size**: Number of ensemble members to generate noise for (must be specified) +- **Temporal resolution** (`dt`): Time step in hours (e.g., `3` for 3-hourly data, `1` for hourly data) +- **Spatial resolution** (`ds`): Grid spacing in km (e.g., `12.5` for 12.5 km grid, `3` for 3 km grid) +- **Spatial sampling** (`nn`): Number of points used for noise generation before upscaling (limited by computational resources) +- **Temporal persistence** (`rho`): AR-1 correlation coefficient (adjust based on forcing time resolution) +- **Covariance matrix**: Cross-variable covariance structure and correlation coefficients (configurable) +- **File paths**: Input/output directories (must be adjusted for local system) + +**Output**: One NetCDF file per month containing noise fields for all ensemble members and all perturbed variables. + +##### Stream File Configuration ##### + +Stream files configure how eCLM reads and applies noise fields to +atmospheric forcing data. Three separate stream files are required +(precipitation, solar radiation, and other variables). Each stream +file must specify the following metadata fields: + +**Required Metadata:** + +- `timeInformation`: Temporal resolution of forcing data in hours (e.g., `3` for 3-hourly data) +- `caseId`: Ensemble member index (e.g., `0` for first member, `1` for second member) +- `numEns`: Total ensemble size used when generating the noise files +- **Perturbation variables**: `tbot_noise` (temperature), `lwdn_noise` (longwave radiation), `precn_noise` (precipitation), `swdn_noise` (shortwave radiation) + + +##### Stream File Examples ##### + +The following examples demonstrate proper stream file configuration +for precipitation, solar radiation, and other atmospheric +variables. Each example must be adapted with correct file paths, +temporal parameters, and ensemble member IDs for your specific +application. + +###### Precipitation Noise Stream File ###### + +**File**: `user_datm.streams.precip_noise.stream_0000.txt` + +```xml + + GENERIC + + + + time time + xc lon + yc lat + area area + mask mask + + + /path/to/domain/file + + + domain.lnd.EUR-11_EUR-11.230216_mask.nc + + + + + PRECTmms precn_noise + + + /path/to/noise/files + + +2003-01.nc +2003-02.nc +2003-03.nc +2003-04.nc +2003-05.nc +2003-06.nc +2003-07.nc +2003-08.nc +2003-09.nc +2003-10.nc +2003-11.nc +2003-12.nc + + + 0 + + + 3 + + + 0 + + + 64 + + +``` + +###### Solar Radiation Noise Stream File ###### + +**File**: `user_datm.streams.solar_noise.stream_0000.txt` + +```xml + + GENERIC + + + + time time + xc lon + yc lat + area area + mask mask + + + /path/to/domain/file + + + domain.lnd.EUR-11_EUR-11.230216_mask.nc + + + + + FSDS swdn_noise + + + /path/to/noise/files + + +2003-01.nc +2003-02.nc +2003-03.nc +2003-04.nc +2003-05.nc +2003-06.nc +2003-07.nc +2003-08.nc +2003-09.nc +2003-10.nc +2003-11.nc +2003-12.nc + + + 0 + + + 3 + + + 0 + + + 64 + + +``` + +###### Temperature and Longwave Radiation Noise Stream File ###### + +**File**: `user_datm.streams.other_noise.stream_0000.txt` + +```xml + + GENERIC + + + + time time + xc lon + yc lat + area area + mask mask + + + /path/to/domain/file + + + domain.lnd.EUR-11_EUR-11.230216_mask.nc + + + + + TBOT tbot_noise + FLDS lwdn_noise + + + /path/to/noise/files + + +2003-01.nc +2003-02.nc +2003-03.nc +2003-04.nc +2003-05.nc +2003-06.nc +2003-07.nc +2003-08.nc +2003-09.nc +2003-10.nc +2003-11.nc +2003-12.nc + + + 0 + + + 3 + + + 0 + + + 64 + + + +``` + +### 3. DATM Integration (`datm_comp_mod.F90`) ### + +Passes ensemble information (`inst_index`, `dt_option`, `ninst`) from +the data atmosphere (DATM) component to the stream infrastructure, +enabling the noise perturbation mechanism to operate correctly within +CESM's multi-instance framework. + +## Design Rationale ## + +1. **Dual perturbation approach**: + - Parameter perturbation for soil properties (offline preprocessing) + - Forcing perturbation via temporal noise patterns (runtime) + +2. **Backward compatibility**: All changes are guarded by `USE_PDAF` + preprocessor directives + +3. **Efficiency**: Uses a single noise file for all ensemble members + by time-shifting, avoiding storage duplication + +4. **Flexibility**: Soil parameters can be perturbed or computed + traditionally based on file contents + +## Modified Source Code ## + +- `src/clm5/biogeophys/SoilStateInitTimeConstMod.F90` +- `src/csm_share/streams/shr_stream_mod.F90` +- `src/csm_share/streams/shr_strdata_mod.F90` +- `src/csm_share/streams/shr_dmodel_mod.F90` +- `src/datm/datm_comp_mod.F90` + diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index a157e16bb6..b3054d47dd 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -150,12 +150,26 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: gti (:) ! read in - fmax real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio) real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio) +#ifdef USE_PDAF + real(r8) ,pointer :: psis_sat (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: shape_param (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: thetas (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: ks (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) + + real(r8) ,pointer :: psis_sat_adj (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: shape_param_adj (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: thetas_adj (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: ks_adj (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) +#endif real(r8) ,pointer :: organic3d (:,:) ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio) character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp integer :: begc, endc integer :: begg, endg +#ifdef USE_PDAF + logical :: parameters_in_file, parameters_in_file_adj +#endif !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp @@ -225,6 +239,17 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) +#ifdef USE_PDAF + allocate(thetas(begg:endg,nlevsoifl)) + allocate(shape_param(begg:endg,nlevsoifl)) + allocate(psis_sat(begg:endg,nlevsoifl)) + allocate(ks(begg:endg,nlevsoifl)) + + allocate(thetas_adj(begg:endg,nlevgrnd)) + allocate(shape_param_adj(begg:endg,nlevgrnd)) + allocate(psis_sat_adj(begg:endg,nlevgrnd)) + allocate(ks_adj(begg:endg,nlevgrnd)) +#endif ! Determine organic_max from parameter file @@ -262,6 +287,60 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if +#ifdef USE_PDAF + ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were + ! generated without parameter perturbation and parameter as input variables + + call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + else + parameters_in_file = .True. + end if + + if (parameters_in_file) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters + ! via pedotransfer function + call ncd_io(ncid=ncid, varname='SHAPE_PARAM', flag='read', data=shape_param, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + + call ncd_io(ncid=ncid, varname='PSIS_SAT', flag='read', data=psis_sat, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + + call ncd_io(ncid=ncid, varname='KSAT', flag='read', data=ks, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + end if + + call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + else + parameters_in_file_adj = .True. + end if + + if (parameters_in_file_adj) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters + ! via pedotransfer function + call ncd_io(ncid=ncid, varname='SHAPE_PARAM_adj', flag='read', data=shape_param_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + + call ncd_io(ncid=ncid, varname='PSIS_SAT_adj', flag='read', data=psis_sat_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + + call ncd_io(ncid=ncid, varname='KSAT_adj', flag='read', data=ks_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + end if +#endif do p = begp,endp g = patch%gridcell(p) if ( sand3d(g,1)+clay3d(g,1) == 0.0_r8 )then @@ -454,6 +533,26 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call pedotransf(ipedof, sand, clay, & soilstate_inst%watsat_col(c,lev), soilstate_inst%bsw_col(c,lev), soilstate_inst%sucsat_col(c,lev), xksat) +#ifdef USE_PDAF + ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer + ! function is used + + if (parameters_in_file) then + if (lev <= nlevsoifl) then + ! Use values from the file for the soil layers + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), lev) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), lev) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), lev) + xksat = ks(col%gridcell(c), lev) ! mm/s + else + ! Use the value from the 10th (default of nlevsoifl) soil level as a default value + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), nlevsoifl) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), nlevsoifl) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), nlevsoifl) + xksat = ks(col%gridcell(c), nlevsoifl) + end if + end if +#endif om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8) om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8) om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8) @@ -462,7 +561,15 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8 soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) +#ifndef USE_PDAF soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b +#else + if (parameters_in_file) then + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b + else + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b + end if +#endif soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat @@ -486,6 +593,15 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat +#ifdef USE_PDAF + if (parameters_in_file_adj) then + ! Use values from the file for the soil layers + soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) + soilstate_inst%bsw_col(c,lev) = shape_param_adj(col%gridcell(c), lev) + soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev) + soilstate_inst%hksat_col(c,lev) = ks_adj(col%gridcell(c), lev) ! mm/s + end if +#endif soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) @@ -559,8 +675,15 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac)*soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) - +#ifndef USE_PDAF soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake +#else + if (parameters_in_file .or. parameters_in_file_adj) then + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake + else + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake + end if +#endif soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac @@ -624,6 +747,10 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) +#ifdef USE_PDAF + deallocate(thetas, shape_param, psis_sat, ks) + deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) +#endif end subroutine SoilStateInitTimeConst diff --git a/src/csm_share/streams/shr_dmodel_mod.F90 b/src/csm_share/streams/shr_dmodel_mod.F90 index 011cb8ccd0..92809378fa 100644 --- a/src/csm_share/streams/shr_dmodel_mod.F90 +++ b/src/csm_share/streams/shr_dmodel_mod.F90 @@ -764,6 +764,13 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs logical :: fileopen character(CL) :: currfile +#ifdef USE_PDAF + ! Yorck: + character(CL) :: mfldName + integer :: position ! integer to look if string contains noise + integer :: caseId ! which ensemble member? + integer :: dt ! time sampling of forcings +#endif integer(in) :: ndims integer(in),pointer :: dimid(:) type(file_desc_t) :: pioid @@ -910,9 +917,42 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs if (my_task == master_task) then call shr_stream_getFileFieldName(stream,k,sfldName) endif +#ifdef USE_PDAF + ! Yorck adjustment + if (my_task == master_task) then + call shr_stream_getModelFieldName(stream,k,mfldName) + end if + call shr_mpi_bcast(mfldName,mpicom,'mfldName') +#endif call shr_mpi_bcast(sfldName,mpicom,'sfldName') rcode = pio_inq_varid(pioid,trim(sfldName),varid) +#ifdef USE_PDAF + !------------------------------------------------------------------------------------ + ! Yorck adjustments : + ! if noise paramter, the frame has to be adjusted depending on the ensemble member + ! it is nt + caseId*24/dt, so that for case 0 (first ensemble member) it takes for + ! the first time step the first frame, for ens mem 1 the 8th (3 hourly data) and + ! so on. With this we ensure different forcings in time and time correlations. + ! idea: plug time and ensemble member into stream + + if (INDEX(mfldName, "noise") == 0) then +#endif frame = nt +#ifdef USE_PDAF + else + if (my_task == master_task) then + call shr_stream_get_dt_and_caseId(stream,dt,caseId) + end if + + call shr_mpi_bcast(caseId,mpicom) + call shr_mpi_bcast(dt,mpicom) + + if (dt == 0) then + call shr_sys_abort(subname//"dt is zero, time resolution of forcing data has to be provided in noise streams") + end if + frame = nt+caseId*24/dt + end if +#endif call pio_setframe(pioid,varid,frame) call pio_read_darray(pioid,varid,pio_iodesc,av%rattr(k,:),rcode) enddo diff --git a/src/csm_share/streams/shr_strdata_mod.F90 b/src/csm_share/streams/shr_strdata_mod.F90 index 3b104d58e7..83c83c8dfc 100644 --- a/src/csm_share/streams/shr_strdata_mod.F90 +++ b/src/csm_share/streams/shr_strdata_mod.F90 @@ -1391,6 +1391,9 @@ subroutine shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, n domTvarName, domXvarName, domYvarName, domAreaName, domMaskName, & filePath, filename, fldListFile, fldListModel, & !--- strdata optional --- +#ifdef USE_PDAF + caseId, dt, numEns, & +#endif nzg, domZvarName, & taxMode, dtlimit, tintalgo, readmode, & fillalgo, fillmask, fillread, fillwrite, & @@ -1441,6 +1444,12 @@ subroutine shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, n character(*),optional ,intent(in) :: readmode ! file read mode character(*),optional, intent(in) :: calendar +#ifdef USE_PDAF + ! Yorck + integer(IN) ,optional ,intent(in) :: caseId + integer(IN) ,optional ,intent(in) :: dt + integer(IN) ,optional ,intent(in) :: numEns +#endif !EOP ! --- local --- @@ -1505,17 +1514,58 @@ subroutine shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, n !---- Backwards compatibility requires Z be optional if (present(domZvarName)) then +#ifndef USE_PDAF call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,taxMode, & fldListFile,fldListModel,domFilePath,domFileName, & domTvarName,domXvarName,domYvarName,domZvarName, & domAreaName,domMaskName, & filePath,filename,trim(name)) +#else + ! Yorck adjustments: optionality of caseId, dt and numEns + if (present(caseId)) then + call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset, & + caseId,dt,numEns, & + taxMode, & + fldListFile,fldListModel,domFilePath,domFileName, & + domTvarName,domXvarName,domYvarName,domZvarName, & + domAreaName,domMaskName, & + filePath,filename,trim(name)) + else + call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset, & + 0,0,0, & + taxMode, & + fldListFile,fldListModel,domFilePath,domFileName, & + domTvarName,domXvarName,domYvarName,domZvarName, & + domAreaName,domMaskName, & + filePath,filename,trim(name)) + end if +#endif else +#ifndef USE_PDAF call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,taxMode, & fldListFile,fldListModel,domFilePath,domFileName, & domTvarName,domXvarName,domYvarName,'undefined', & domAreaName,domMaskName, & filePath,filename,trim(name)) +#else + if (present(caseId)) then + call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset, & + caseId,dt,numEns, & + taxMode, & + fldListFile,fldListModel,domFilePath,domFileName, & + domTvarName,domXvarName,domYvarName,'undefined', & + domAreaName,domMaskName, & + filePath,filename,trim(name)) + else + call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset, & + 0,0,0, & + taxMode, & + fldListFile,fldListModel,domFilePath,domFileName, & + domTvarName,domXvarName,domYvarName,'undefined', & + domAreaName,domMaskName, & + filePath,filename,trim(name)) + end if +#endif endif if (present(nzg)) then @@ -1536,6 +1586,9 @@ subroutine shr_strdata_create_oldway(SDAT, name, mpicom, compid, gsmap, ggrid, n filePath, filename, fldListFile, fldListModel, & pio_subsystem, pio_iotype, & !--- strdata optional --- +#ifdef USE_PDAF + caseId, dt, numEns, & +#endif nzg, domZvarName, & taxMode, dtlimit, tintalgo, readmode, & fillalgo, fillmask, fillread, fillwrite, & @@ -1589,8 +1642,15 @@ subroutine shr_strdata_create_oldway(SDAT, name, mpicom, compid, gsmap, ggrid, n character(*),optional ,intent(in) :: readmode ! file read mode character(*),optional, intent(in) :: calendar +#ifdef USE_PDAF + ! Yorck + integer(IN) ,optional ,intent(in) :: caseId + integer(IN) ,optional ,intent(in) :: dt + integer(IN) ,optional ,intent(in) :: numEns +#endif !EOP ! pio variables are already in SDAT no need to copy them +#ifndef USE_PDAF call shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, nxg, nyg,& yearFirst, yearLast, yearAlign, offset, domFilePath, domFilename, domTvarName, & domXvarName, domYvarName, domAreaName, domMaskName, & @@ -1598,7 +1658,17 @@ subroutine shr_strdata_create_oldway(SDAT, name, mpicom, compid, gsmap, ggrid, n taxMode=taxMode,dtlimit=dtlimit,tintalgo=tintalgo,readmode=readmode,& fillalgo=fillalgo,fillmask=fillmask,fillread=fillread,fillwrite=fillwrite,mapalgo=mapalgo,& mapmask=mapmask,mapread=mapread,mapwrite=mapwrite,calendar=calendar) - +#else + call shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, nxg, nyg,& + yearFirst, yearLast, yearAlign, offset, domFilePath, domFilename, domTvarName, & + domXvarName, domYvarName, domAreaName, domMaskName, & + filePath, filename, fldListFile, fldListModel, & + caseId=caseId, dt=dt, numEns=numEns, & + nzg=nzg, domZvarName=domZvarName,& + taxMode=taxMode,dtlimit=dtlimit,tintalgo=tintalgo,readmode=readmode,& + fillalgo=fillalgo,fillmask=fillmask,fillread=fillread,fillwrite=fillwrite,mapalgo=mapalgo,& + mapmask=mapmask,mapread=mapread,mapwrite=mapwrite,calendar=calendar) +#endif end subroutine shr_strdata_create_oldway !=============================================================================== diff --git a/src/csm_share/streams/shr_stream_mod.F90 b/src/csm_share/streams/shr_stream_mod.F90 index 1e8d5c0bed..4fc371f3a9 100644 --- a/src/csm_share/streams/shr_stream_mod.F90 +++ b/src/csm_share/streams/shr_stream_mod.F90 @@ -77,6 +77,10 @@ module shr_stream_mod public :: shr_stream_setAbort ! set internal shr_stream abort flag public :: shr_stream_getDebug ! get internal shr_stream debug level public :: shr_stream_isInit ! check if stream is initialized +#ifdef USE_PDAF + ! Yorck + public :: shr_stream_get_dt_and_caseId ! return dt and caseId +#endif ! public :: shr_stream_bcast ! broadcast a stream (untested) ! !PUBLIC DATA MEMBERS: @@ -151,6 +155,16 @@ module shr_stream_mod character(SHR_KIND_CS) :: domAreaName ! domain file: area var name character(SHR_KIND_CS) :: domMaskName ! domain file: mask var name +#ifdef USE_PDAF + ! Yorck: in the stream, also the time resolution of the forcing data and the caseId + ! is saved. This is done to get the right frame of the noise file later. + ! Also: the number of ensemble members the noise data was produced for has + ! to be included to make sure that the time information is right + integer(SHR_KIND_IN) :: caseId + integer(SHR_KIND_IN) :: dt + integer(SHR_KIND_IN) :: numEns +#endif + character(SHR_KIND_CS) :: tInterpAlgo ! Algorithm to use for time interpolation character(SHR_KIND_CL) :: calendar ! stream calendar end type shr_stream_streamType @@ -642,6 +656,78 @@ subroutine shr_stream_init(strm,infoFile,yearFirst,yearLast,yearAlign,taxMode,rc close(nUnit) +#ifdef USE_PDAF + !----------------------------------------------------------------------------- + ! Yorck adaptations: read in time resolution of forcings and caseId, numEns + !----------------------------------------------------------------------------- + + !--- find start tag --- + open(nUnit,file=infoFile,STATUS='OLD',FORM='FORMATTED',ACTION='READ') + startTag = "" + endTag = "" + call shr_stream_readUpToTag(nUnit,startTag,rc=rCode2) + if (rCode2 /= 0)then + rCode = rCode2 + goto 999 + end if + startTag = "" + endTag = "" + call shr_stream_readUpToTag(nUnit,startTag,optionalTag=.true.,rc=rCode2) + if (rCode2 == 0) then + !--- read data --- + read(nUnit,*,END=999) int + strm%caseId = int + else + strm%caseId = 0 + end if + + close(nUnit) + + !--- find start tag --- + open(nUnit,file=infoFile,STATUS='OLD',FORM='FORMATTED',ACTION='READ') + startTag = "" + endTag = "" + call shr_stream_readUpToTag(nUnit,startTag,rc=rCode2) + if (rCode2 /= 0)then + rCode = rCode2 + goto 999 + end if + startTag = "" + endTag = "" + call shr_stream_readUpToTag(nUnit,startTag,optionalTag=.true.,rc=rCode2) + if (rCode2 == 0) then + !--- read data --- + read(nUnit,*,END=999) int + strm%dt = int + else + strm%dt = 0 + end if + + close(nUnit) + + !--- find start tag --- + open(nUnit,file=infoFile,STATUS='OLD',FORM='FORMATTED',ACTION='READ') + startTag = "" + endTag = "" + call shr_stream_readUpToTag(nUnit,startTag,rc=rCode2) + if (rCode2 /= 0)then + rCode = rCode2 + goto 999 + end if + startTag = "" + endTag = "" + call shr_stream_readUpToTag(nUnit,startTag,optionalTag=.true.,rc=rCode2) + if (rCode2 == 0) then + !--- read data --- + read(nUnit,*,END=999) int + strm%numEns = int + else + strm%numEns = 0 + end if + + close(nUnit) +#endif + !----------------------------------------------------------------------------- ! get initial calendar value !----------------------------------------------------------------------------- @@ -681,12 +767,19 @@ end subroutine shr_stream_init ! ! !INTERFACE: ------------------------------------------------------------------ +#ifdef USE_PDAF + subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset, caseId, dt, & + numEns, taxMode, fldListFile,fldListModel,domFilePath,domFileName, & + domTvarName,domXvarName,domYvarName,domZvarName, & + domAreaName,domMaskName, & + filePath,filename,dataSource,rc) +#else subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, & fldListFile,fldListModel,domFilePath,domFileName, & domTvarName,domXvarName,domYvarName,domZvarName, & domAreaName,domMaskName, & filePath,filename,dataSource,rc) - +#endif ! !INPUT/OUTPUT PARAMETERS: type(shr_stream_streamType) ,intent(inout) :: strm ! data stream @@ -710,6 +803,12 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, & character(*) ,optional,intent(in) :: dataSource ! comment line integer (SHR_KIND_IN),optional,intent(out) :: rc ! return code +#ifdef USE_PDAF + ! Yorck + integer (SHR_KIND_IN),optional,intent(in) :: caseId ! ensemble member case ID + integer (SHR_KIND_IN),optional,intent(in) :: dt ! time resolution of forcing data + integer (SHR_KIND_IN),optional,intent(in) :: numEns ! number of ensemble members set to produce noise data +#endif !EOP !----- local ----- @@ -791,6 +890,21 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, & call fileVec%move_out(strm%file) endif +#ifdef USE_PDAF + ! Yorck + if (present(caseId)) then + strm%caseId = caseId + end if + + if (present(dt)) then + strm%dt = dt + end if + + if (present(numEns)) then + strm%numEns = numEns + end if +#endif + !----------------------------------------------------------------------------- ! get initial calendar value !----------------------------------------------------------------------------- @@ -871,6 +985,12 @@ subroutine shr_stream_default(strm,rc) strm%domAreaName = ' ' strm%domMaskName = ' ' +#ifdef USE_PDAF + strm%caseId = 0 + strm%dt = 0 + strm%numEns = 0 +#endif + strm%calendar = shr_cal_noleap if ( present(rc) ) rc = 0 @@ -1516,6 +1636,11 @@ subroutine shr_stream_readTCoord(strm,k,rc) integer(SHR_KIND_IN) :: ndate ! calendar date of time value real(SHR_KIND_R8) :: nsec ! elapsed secs on calendar date real(SHR_KIND_R8),allocatable :: tvar(:) +#ifdef USE_PDAF + character(SHR_KIND_CL) :: mfldName !Yorck + integer(SHR_KIND_IN) :: nt_ ! Yorck + real(SHR_KIND_R8),allocatable :: tvar_lok(:) ! Yorck +#endif !----- formats ----- character(*),parameter :: subname = '(shr_stream_readTCoord) ' character(*),parameter :: F01 = "('(shr_stream_readTCoord) ',a,2i7)" @@ -1539,6 +1664,28 @@ subroutine shr_stream_readTCoord(strm,k,rc) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension') deallocate(dids) +#ifdef USE_PDAF + ! Yorck + ! as the time information in the noise files is wrong, the last timesteps which only resemble the number of ensemble + ! members and the time resolution of the forcing data have to be deleted. It is not necessarily the number + ! of ensemble members the experiment is conducted with but the number of ensemble members the noise data was produced with. + + ! After nt is adjusted, also the read in number of timesteps has to be adjusted. + + call shr_stream_getModelFieldName(strm,1,mfldName) + if (INDEX(mfldName, "noise") .ne. 0) then + nt_ = nt + if (strm%numEns == 0) then + call shr_sys_abort(subname//' ERROR: number of ensemble members that were used for producing noise data has to be set in stream') + end if + if (strm%dt == 0) then + call shr_sys_abort(subname//' ERROR: time resolution of forcing data has to be set in stream') + end if + ! (numEns-1)*(24/dt) is the number of extra time steps introduced in the noise preperation + nt = nt - (strm%numEns-1)*(24/strm%dt) + end if +#endif + allocate(strm%file(k)%date(nt),strm%file(k)%secs(nt)) strm%file(k)%nt = nt @@ -1563,10 +1710,30 @@ subroutine shr_stream_readTCoord(strm,k,rc) strm%calendar = trim(shr_cal_calendarName(trim(calendar))) allocate(tvar(nt)) +#ifdef USE_PDAF + ! Yorck + if (INDEX(mfldName, "noise") .ne. 0) then + allocate(tvar_lok(nt_)) + rcode = nf90_get_var(fid,vid,tvar_lok) + if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var') + rCode = nf90_close(fid) + if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_close') + + do n = 1, nt + tvar(n) = tvar_lok(n) + end do + deallocate(tvar_lok) + + else +#endif rcode = nf90_get_var(fid,vid,tvar) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var') rCode = nf90_close(fid) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_close') +#ifdef USE_PDAF + end if + ! End Yorck +#endif do n = 1,nt call shr_cal_advDate(tvar(n),bunits,bdate,bsec,ndate,nsec,calendar) strm%file(k)%date(n) = ndate @@ -2582,6 +2749,11 @@ subroutine shr_stream_restWrite(strm,fileName,caseName,caseDesc,nstrms,rc) write(nUnit) strm(k)%domAreaName ! domain file: area var name write(nUnit) strm(k)%domMaskName ! domain file: mask var name +#ifdef USE_PDAF + write(nUnit) strm(k)%dt ! Yorck + write(nUnit) strm(k)%caseId ! Yorck + write(nUnit) strm(k)%numEns ! Yorck +#endif end do close(nUnit) @@ -2885,6 +3057,12 @@ subroutine shr_stream_dataDump(strm) write(s_logunit,F00) "domAreaName = ", trim(strm%domAreaName) write(s_logunit,F00) "domMaskName = ", trim(strm%domMaskName) +#ifdef USE_PDAF + write(s_logunit,F01) "dt = ", strm%dt ! Yorck + write(s_logunit,F01) "caseId = ", strm%caseId ! Yorck + write(s_logunit,F01) "numEns = ", strm%numEns !Yorck +#endif + end subroutine shr_stream_dataDump !=============================================================================== @@ -3232,6 +3410,12 @@ subroutine shr_stream_bcast(stream,comm,rc) call shr_mpi_bcast(stream%domZvarName ,comm,subName) call shr_mpi_bcast(stream%domMaskName ,comm,subName) call shr_mpi_bcast(stream%calendar ,comm,subName) + +#ifdef USE_PDAF + !Yorck + call shr_mpi_bcast(stream%dt ,comm,subName) + call shr_mpi_bcast(stream%caseId ,comm,subName) +#endif if (pid /= 0) allocate(stream%file(stream%nFiles)) @@ -3247,6 +3431,31 @@ subroutine shr_stream_bcast(stream,comm,rc) end subroutine shr_stream_bcast +#ifdef USE_PDAF + ! Yorck + subroutine shr_stream_get_dt_and_caseId(stream,dt,caseId) + + implicit none + + ! !INPUT/OUTPUT PARAMETERS: + + type(shr_stream_streamType) ,intent(in) :: stream ! stream in question + integer(SHR_KIND_IN) ,intent(out) :: dt ! + integer(SHR_KIND_IN) ,intent(out) :: caseId ! + + + + !------------------------------------------------------------------------------- + ! Notes: + !------------------------------------------------------------------------------- + + dt = stream%dt + caseId = stream%caseId + + + end subroutine shr_stream_get_dt_and_caseId +#endif + !=============================================================================== end module shr_stream_mod !=============================================================================== diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index 0b7744d1da..25f4b6b018 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/datm/datm_comp_mod.F90 @@ -77,6 +77,12 @@ module datm_comp_mod integer(IN) :: krc,krl,ksc,ksl,kswndr,kswndf,kswvdr,kswvdf,kswnet integer(IN) :: kanidr,kanidf,kavsdr,kavsdf integer(IN) :: stbot,swind,sz,spbot,sshum,stdew,srh,slwdn,sswdn,sswdndf,sswdndr +#ifdef USE_PDAF + integer(IN) :: stbot_noise ! temperature noise + integer(IN) :: sprecn_noise ! precipitation noise + integer(IN) :: sswdn_noise ! shortwave radiation noise + integer(IN) :: slwdn_noise ! longwave radiation noise +#endif integer(IN) :: sprecc,sprecl,sprecn,sco2p,sco2d,sswup,sprec,starcf ! water isotopes / tracer input @@ -161,7 +167,11 @@ module datm_comp_mod ! other fields used in calculations. Fields that are simply read and passed directly to ! the coupler do not need to be in these lists. +#ifdef USE_PDAF + integer(IN),parameter :: ktranss = 37 +#else integer(IN),parameter :: ktranss = 33 +#endif character(16),parameter :: stofld(1:ktranss) = & (/"strm_tbot ","strm_wind ","strm_z ","strm_pbot ", & @@ -175,6 +185,10 @@ module datm_comp_mod "strm_pbot_af ","strm_shum_af ","strm_swdn_af ","strm_lwdn_af ", & "strm_rh_18O ","strm_rh_HDO ", & "strm_precn_16O ","strm_precn_18O ","strm_precn_HDO " & +#ifdef USE_PDAF + ! add DATM-internal stream variable names for noise fields (Yorck) + ,"strm_tbot_noise ","strm_precn_noise","strm_swdn_noise ","strm_lwdn_noise " & +#endif /) character(16),parameter :: stifld(1:ktranss) = & @@ -190,6 +204,10 @@ module datm_comp_mod ! isotopic forcing "rh_18O ","rh_HDO ", & "precn_16O ","precn_18O ","precn_HDO " & +#ifdef USE_PDAF + ! add NetCDF input file variable names for noise fields (Yorck) + ,"tbot_noise ","precn_noise ","swdn_noise ","lwdn_noise " & +#endif /) character(CL), pointer :: ilist_av(:) ! input list for translation @@ -444,6 +462,14 @@ subroutine datm_comp_init(Eclock, x2a, a2x, & sprec = mct_aVect_indexRA(avstrm,'strm_prec' ,perrWith='quiet') starcf = mct_aVect_indexRA(avstrm,'strm_tarcf' ,perrWith='quiet') +#ifdef USE_PDAF + ! Get Noise Field Indices from Stream (noise streams are optional) + stbot_noise = mct_aVect_indexRA(avstrm,'strm_tbot_noise' ,perrWith='quiet') + sprecn_noise = mct_aVect_indexRA(avstrm,'strm_precn_noise' ,perrWith='quiet') + sswdn_noise = mct_aVect_indexRA(avstrm,'strm_swdn_noise' ,perrWith='quiet') + slwdn_noise = mct_aVect_indexRA(avstrm,'strm_lwdn_noise' ,perrWith='quiet') +#endif + ! anomaly forcing sprecsf = mct_aVect_indexRA(avstrm,'strm_precsf' ,perrWith='quiet') sprec_af = mct_aVect_indexRA(avstrm,'strm_prec_af' ,perrWith='quiet') @@ -902,6 +928,14 @@ subroutine datm_comp_run(EClock, x2a, a2x, & endif if (my_task == master_task) & write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax +#ifdef USE_PDAF + ! Log if temperature noise read from stream file + if (my_task == master_task) then + if (stbot_noise > 1 ) then + write(logunit,*) trim(subname),' Using noise from files ' + end if + end if +#endif endif lsize = mct_avect_lsize(a2x) do n = 1,lsize @@ -910,6 +944,40 @@ subroutine datm_comp_run(EClock, x2a, a2x, & !--- temperature --- if (tbotmax < 50.0_R8) a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + tkFrz +#ifdef USE_PDAF + !!! begin perturbation block (Yorck) --> perturb + !!! temperature, long wave radiation, short wave radiation + !!! and precipitation + ! --> further variables can be introduced (changes also have + ! to be made in reading in routines of streams for this, for + ! now only these variables are included) + + ! tbot + if (stbot_noise > 0) then + ! Additive noise + a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + avstrm%rAttr(stbot_noise,n) + end if + + ! lwdn + if (slwdn_noise > 0) then + ! Additive noise + avstrm%rAttr(slwdn,n) = avstrm%rAttr(slwdn,n) + avstrm%rAttr(slwdn_noise,n) + end if + + ! swdn + if (sswdn_noise>0) then + ! Multiplicative noise + avstrm%rAttr(sswdn,n) = avstrm%rAttr(sswdn,n) * avstrm%rAttr(sswdn_noise,n) + end if + + ! sprecn + if (sprecn_noise > 0) then + ! Multiplicative noise + avstrm%rAttr(sprecn,n) = avstrm%rAttr(sprecn,n)*avstrm%rAttr(sprecn_noise,n) + end if + + !!! end perturbation block +#endif ! Limit very cold forcing to 180K a2x%rAttr(ktbot,n) = max(180._r8, a2x%rAttr(ktbot,n)) a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)