diff --git a/fer/gnl/load_dsg_mask_ftrset_var.F b/fer/gnl/load_dsg_mask_ftrset_var.F new file mode 100644 index 000000000..4442ac343 --- /dev/null +++ b/fer/gnl/load_dsg_mask_ftrset_var.F @@ -0,0 +1,225 @@ + SUBROUTINE LOAD_DSG_MASK_FTRSET_VAR (dset, mask_var, reload, status) + +* set data/trajmask= or set data/stnmask= +* +* If /SMASK is used to define a trajectory-mask was given for trajectoryProfile data, or +* a timeseries-station mask for timeseriesProfile data, store the values of the mask in +* line storage.Save the memory location dsg_ftrsetmsk_lm(dset), in common. Mask data used +* in creating feature-masks when computing feature-level masking which also may include +* context limits in the E direftion. Will deallocate line mem dsg_ftrsetmsk_lm on +* defining a new mask, closing the dataset, or CANCEL DATA/FMASK. +* If there is a /FMASK then cancel it; if there is a /SMASK then +* it will be canceled upon defining an /FMASK. +* +* On a "LET mask= " find dataset(s) where the mask being defined has +* the same name as the dataset's current mask. Here also test the +* length of the new mask. If it's the right length for this dataset, +* will replace its FMASK, otherwise the FMASK is canceled. +* +* 8/20/2020 *acm* + + include 'tmap_dims.parm' + include 'ferret.parm' + include 'errmsg.parm' + include 'slash.parm' + include 'command.parm' + include 'xvariables.cmn' + include 'xcontext.cmn' + include 'xprog_state.cmn' +#include "tmap_dset.parm" + include 'xdset_info.cmn_text' + include 'xtm_grid.cmn_text' + include 'netcdf.inc' + + LOGICAL reload + INTEGER dset, status + CHARACTER*(*) mask_var + + LOGICAL TM_LEGAL_NAME, TM_ITSA_DSG, is_cmpnd + INTEGER TM_DSG_DSET_FROM_GRID, TM_DSG_NF2FEATURES, + . TM_LENSTR1, CX_DIM_LEN, + . i, ivar, grid, loc, ndim, idir, npts, cx, mr1, + . nftrsets, idim, dim(nferdims), attlen, varid, + . attid, attype, attoutflag, slen, orient + CHARACTER FULL_VAR_TITLE*128, TM_FMT*48, + . buff*128, mask_title*128, buff2*7 + REAL dummy + +c Make sure the name is ok + + IF (mask_var .EQ. ' ') GOTO 5100 + IF ( .NOT.TM_LEGAL_NAME(mask_var) ) GOTO 5200 + +* For error messages and notes + buff2 = 'station' + IF (orient .EQ. pfeatureType_TrajectoryProfile) buff2 = 'traject' + +* Is it a timeseriesProfile or trajectoryProfile dataset? + + orient = dsg_orientation(dset) + is_cmpnd = orient .EQ. pfeatureType_TrajectoryProfile .OR. + . orient .EQ. pfeatureType_TimeseriesProfile + + IF (.NOT.is_cmpnd) GOTO 5300 + + nftrsets = TM_DSG_NF2FEATURES (dset) + +* Set the line-memory to store the mask data + +* ... A feature-mask was already applied to the dataset. Wipe that one out. + + IF (dsg_msk_lm(dset) .NE. unspecified_int4 .OR. + . dsg_ftrsetmsk_lm(dset) .NE. unspecified_int4) THEN + + IF (dsg_ftrsetmsk_lm(dset) .NE. unspecified_int4) THEN + CALL FREE_LINE_DYNMEM( dsg_ftrsetmsk_lm(dset) ) + CALL TM_DEALLO_DYN_LINE( dsg_ftrsetmsk_lm(dset) ) + dsg_ftrsetmsk_lm(dset) = unspecified_int4 + ENDIF + IF (dsg_msk_lm(dset) .NE. unspecified_int4) THEN + CALL FREE_LINE_DYNMEM( dsg_msk_lm(dset) ) + CALL TM_DEALLO_DYN_LINE( dsg_msk_lm(dset) ) + dsg_msk_lm(dset) = unspecified_int4 + ENDIF + +* ... wipe memory clear of stored variables - this could change all definitions + DO i = 1,max_mr_avail + IF ( mr_protected( i ) .NE. mr_deleted ) + . CALL DELETE_VARIABLE( i ) + ENDDO + + ENDIF + +* Load the variable. Reset the arg pointers to get the mask-var +* Use arg_start, arg_end to point to the name of the variable. +* We may be coming from a SET DATA/TRAJMASK= or SET DATA/STNMASK= +* command or a LET command redefining a variable whose name is +* already associated with a feature-mask for this dataset. + + loc = qual_given(slash_set_data_smask) + IF (loc.GT.0) THEN + IF (.NOT. is_cmpnd) GOTO 5400 + ENDIF + + IF (loc .GT. 0) THEN + num_args = 1 + arg_end(1) = qual_end(loc) + arg_start(1) = qual_end(loc) - TM_LENSTR1(mask_var) + 1 + ELSE + arg_end(1) = arg_start(1) + TM_LENSTR1(mask_var) - 1 + ENDIF + + CALL GET_CMND_DATA ( cx_last, ptype_float, status ) + IF (status .NE. ferr_ok .OR. num_uvars_in_cmnd .GT. 1) GOTO 5500 + +* ... make sure it's a line (not a point,plane,etc.) and that the length +* is the nftrsets-length: the # of stations or # of trajectories + + cx = is_cx(isp) + CALL GET_CX_DIMS( cx, ndim, dim ) + IF ( ndim .GT. 1 ) GOTO 5500 + idim = dim(1) + IF (nftrsets .NE. CX_DIM_LEN( idim, cx ) ) THEN +* the mask doesn't match mask on this dset. Issue a note + IF (reload) GOTO 5600 + GOTO 5500 + ENDIF + +* store it here, + + CALL TM_ALLO_TMP_LINE(dsg_ftrsetmsk_lm(dset), status) + CALL GET_LINE_DYNMEM (nftrsets, dsg_ftrsetmsk_lm(dset), status) + CALL TM_USE_LINE(dsg_ftrsetmsk_lm(dset)) + + CALL TM_NEW_LINE_NAME ('FEATURE_MASK', buff ) + line_name(dsg_ftrsetmsk_lm(dset)) = buff + line_direction( dsg_ftrsetmsk_lm(dset) ) = 'EE' + + mr1 = is_mr( num_uvars_in_cmnd ) + cx = is_cx( num_uvars_in_cmnd ) + + CALL EXTRACT_LINE (cx, + . memry(mr1)%ptr, + . mr1, + . linemem(dsg_ftrsetmsk_lm(dset))%ptr, + . idir, + . ndim, + . npts, + . status ) + + IF (npts.NE.nftrsets .OR. status.NE.ferr_ok) GOTO 5000 + +* Put the default bad-flag into the mask data + CALL TM_SWITCH_BAD ( cx_bad_data (cx), bad_val4, linemem(dsg_ftrsetmsk_lm(dset))%ptr, npts ) + +* Set a global attribute with the mask title + + buff = '__feature_mask_' + attlen = TM_LENSTR1(buff) + mask_title = FULL_VAR_TITLE( cx, .FALSE., attlen ) + + varid = 0 + CALL CD_GET_VAR_ATT_ID (dset, varid, buff, attid, status) + IF (attid .GT. 0) THEN + CALL CD_GET_VAR_ATT_INFO (dset, varid, attid, + . buff, attype, attlen, attoutflag, status ) + attoutflag = 0 + CALL CD_REPLACE_ATTR (dset, varid, buff, attype, attlen, + . mask_title, dummy, status) + ELSE + attype = NCCHAR + attoutflag = 0 + CALL CD_PUT_NEW_ATTR (dset, varid, buff, attype, + . attlen, attoutflag, mask_title, dummy, status ) + ENDIF + IF (status .NE. ferr_ok) GOTO 5000 + +* And another with the feature-mask name + + buff = '__feature_mask_var' + attlen = TM_LENSTR1(buff) + + varid = 0 + CALL CD_GET_VAR_ATT_ID (dset, varid, buff, attid, status) + IF (attid .GT. 0) THEN + CALL CD_GET_VAR_ATT_INFO (dset, varid, attid, + . buff, attype, attlen, attoutflag, status ) + attoutflag = 0 + CALL CD_REPLACE_ATTR (dset, varid, buff, attype, attlen, + . mask_var, dummy, status) + ELSE + attype = NCCHAR + attoutflag = 0 + CALL CD_PUT_NEW_ATTR (dset, varid, buff, attype, + . attlen, attoutflag, mask_var, dummy, status ) + ENDIF + IF (status .NE. ferr_ok) GOTO 5000 + + status = ferr_ok + + 5000 RETURN + + 5100 CALL ERRMSG( ferr_unknown_arg, status, 'argument required /FMASK=?', *5000) + 5200 CALL ERRMSG( ferr_invalid_command, status, + . '/MASK=name is not an acceptable name', *5000) + 5300 CALL ERRMSG (ferr_invalid_command, status, + .'/'//buff2//' is set only for Discrete Sampling Geometries datasets', + . *5000 ) + 5400 CALL ERRMSG (ferr_invalid_command, status, + .'/SMASK is set only for Discrete Sampling Geometries '//pCR// + .'or TrajectoryProfile datasets', *5000 ) + + 5500 buff = TM_FMT(DBLE(nftrsets), 0, 12, i) + CALL ERRMSG (ferr_invalid_command, status, + .'Mask variable must be 1-dimensional, with length num-'//buff2//' = '// + . buff(:i), *5000 ) + + 5600 buff = TM_FMT(DBLE(dset), 0, 12, i) + slen = TM_LENSTR1(mask_var) + CALL WARN (buff2//'-mask on dataset '//buff(:i)// + . ' canceled. New definition of '//mask_var(:slen)// + . ' length does not match # of '//buff2//' in dataset') + GOTO 5000 + + END + diff --git a/fer/gnl/reload_dsg_ftrmaskvar.F b/fer/gnl/reload_dsg_ftrmaskvar.F new file mode 100644 index 000000000..ecb510d39 --- /dev/null +++ b/fer/gnl/reload_dsg_ftrmaskvar.F @@ -0,0 +1,74 @@ + SUBROUTINE RELOAD_DSG_FTRMASKVAR (varname, slen) + + IMPLICIT NONE +# include "tmap_dset.parm" + include 'tmap_dims.parm' + include 'ferret.parm' + include 'errmsg.parm' + include 'xdset_info.cmn_text' + include 'xprog_state.cmn' + +* 6/3/2020 *acm* +* If a variable that is a station-mask or trajectory-mask on a timeseriesProfile or +* trajectoryProfile dataset is redefined, the previous feature mask info is cleared out. +* Here, check the variable and load the new definition as the feature mask. + + INTEGER slen + CHARACTER*(*) varname + + LOGICAL NC_GET_ATTRIB, got_it, reload + INTEGER STR_SAME, iset, varid, attid, attlen, + . attoutflag, maxlen, alen, status + REAL attval + CHARACTER*32 attname, attstring, dummy + CHARACTER*48 TM_FMT + +* ... loop over datasets. See if the variable is listed as a feature-mask +* for the dataset. If so re-load the mask. + + maxlen = 32 + varid = 0 + + DO iset = pdset_irrelevant+1, maxdsets + + IF ( ds_name(iset) .EQ. char_init2048) CYCLE + + varid = 0 ! look in global attributes for '__feature_mask_var' + attname = '__feature_mask_var' + + CALL CD_GET_VAR_ATT_ID (iset, varid, attname, attid, status) + + IF (status .NE. ferr_ok) CYCLE + got_it = NC_GET_ATTRIB ( iset, varid, attname,.FALSE., dummy, + . maxlen, attlen, attoutflag, attstring, + . attval ) + + IF (.NOT. got_it) CYCLE + IF (STR_SAME(varname(:slen), attstring(:attlen)) .EQ. 0) THEN + +* The mask has the same name as this dataset's current smask. The LOAD +* routine also tests the length of the new mask. If it's the right +* length for this dataset, will replace its SMASK, otherwise the +* SMASK is canceled. + +* These attributes will get restored using the new definition for +* the mask variable in the load routine + + CALL CD_DELETE_ATTRIBUTE (iset, varid, attname, status ) + attname = '__feature_mask_' + CALL CD_GET_VAR_ATT_ID (iset, varid, attname, attid, status) + CALL CD_DELETE_ATTRIBUTE (iset, varid, attname, status ) + + cmnd_buff = 'load '//varname(:slen) + + arg_start(1) = 6 + arg_end(1) = 5+slen + + reload = .TRUE. + CALL LOAD_DSG_MASK_FTRSET_VAR (iset, varname, reload, status) + ENDIF + + ENDDO + + RETURN + END diff --git a/fer/gnl/repl_exprns.F b/fer/gnl/repl_exprns.F index c9113726c..98c11a4d6 100644 --- a/fer/gnl/repl_exprns.F +++ b/fer/gnl/repl_exprns.F @@ -133,6 +133,7 @@ SUBROUTINE REPL_EXPRNS( cmnd, lencmnd, cmnd_num, * correct info for `var,return=dset` and dsetnum and dsetpath. * V745 4/19 *acm* ticket 1916, special handling getting context for RETURN= outputs * V751 6/19 *acm* ticket 1929, don't chop off str_len when testing for reformatting at the end +* V7.63 9/20 *acm* add `var,RETURN=*coord` to get coordinate variables for DSG datasets IMPLICIT NONE #include "netcdf.inc" diff --git a/fer/plt/dsg_pltalong_setup.F b/fer/plt/dsg_pltalong_setup.F new file mode 100644 index 000000000..6e0343195 --- /dev/null +++ b/fer/plt/dsg_pltalong_setup.F @@ -0,0 +1,166 @@ + SUBROUTINE DSG_PLTALONG_SETUP (dset, plot_orient, its_traj, + . dsg_as_traj, dsg_as_time, status) + +* Get PLOT/ALONG= and return the orientation to use in plotting for a DSG dataset +* +* TrajectoryProfile PLOT/ALONG=xy ==> trajectory (dsg_as_traj is set) +* PLOT/ALONG=z ==> prof +* TimeseriesProfile PLOT/ALONG=xy ==> point +* PLOT/ALONG=z ==> prof +* PLOT/ALONG=t ==> time (dsg_as_time is set) +* Profile PLOT/ALONG=xy ==> point +* Timeseries PLOT/ALONG=xy ==> point +* Trajectory PLOT/ALONG=t ==> time +* +* as well, all of the types can have along= to their native +* direction, e.g. Timeseries can have /along=t, etc +* + + IMPLICIT NONE + include 'ferret.parm' + include 'errmsg.parm' + include 'slash.parm' + include 'tmap_dims.parm' + include 'xdset_info.cmn_text' + include 'xprog_state.cmn' + include 'xtext_info.cmn' + +#include "tmap_dset.parm" + +* arguments9/9/2020 + INTEGER dset, plot_orient, status + LOGICAL its_traj, dsg_as_traj, dsg_as_time + +* Local variables + INTEGER orientation, loc, along_dim + CHARACTER buff1*1, message*128 + +* Initialize + message = 'none' + status = ferr_ok + + dsg_as_time = .FALSE. + dsg_as_traj = .FALSE. + + IF (dset .LE. pdset_irrelevant) RETURN + + its_traj = dsg_feature_type(dset).EQ.pfeatureType_Trajectory .OR. + . dsg_feature_type(dset).EQ.pfeatureType_Point + +* dataset orientation. + + orientation = dsg_orientation(dset) + +* Is there a PLOT/ALONG= qualifier? + + along_dim = no_dim + loc = qual_given(slash_plot_along) + IF ( loc .GT. 0. ) THEN + + CALL EQUAL_STRING( cmnd_buff(qual_start(loc):qual_end(loc)), + . buff1, status ) + IF (status .NE. ferr_ok) THEN + plot_orient = orientation + GOTO 4000 + ENDIF + + DO along_dim = 1, nferdims + IF ( buff1 .EQ. ww_dim_name(along_dim) ) EXIT + ENDDO + + ENDIF + + IF (along_dim .EQ. no_dim) GOTO 4000 + +* Check valid setings for PLOT/ALONG= for each Feature Type. +* When the setting agrees with the default plot direction, ignore +* the /along= setting. + + IF (orientation .EQ. pfeatureType_Point) THEN + + plot_orient = orientation + IF (along_dim .EQ. x_dim) THEN + GOTO 5000 + ELSE + message = 'PLOT/ALONG='//ww_dim_name(along_dim)// + . ' : Point data can only be plotted along XY' + ENDIF + ELSE IF (orientation .EQ. pfeatureType_Trajectory) THEN + + IF (along_dim .EQ. x_dim) THEN + GOTO 5000 + ELSE IF (along_dim.EQ.t_dim) THEN + plot_orient = pfeatureType_TimeSeries + ELSE + message = 'PLOT/ALONG='//ww_dim_name(along_dim)// + . ' : Trajectory data can only be plotted along XY or T' + ENDIF + + + ELSE IF (orientation .EQ. pfeatureType_Profile) THEN + + IF (along_dim .EQ. z_dim) THEN + GOTO 5000 + ELSE IF (along_dim .EQ. x_dim) THEN + plot_orient = pfeatureType_Point + dsg_as_traj = .TRUE. + ELSE + message = 'PLOT/ALONG='//ww_dim_name(along_dim)// + . ' : Profile data can only be plotted along XY or Z' + ENDIF + + ELSE IF (orientation .EQ. pfeatureType_TimeSeries) THEN + + IF (along_dim .EQ. t_dim) THEN + GOTO 5000 + ELSE IF (along_dim .EQ. x_dim) THEN + plot_orient = pfeatureType_Point + dsg_as_traj = .TRUE. + ELSE + message = 'PLOT/ALONG='//ww_dim_name(along_dim)// + . ' : Timeseries data can only be plotted along XY or Z' + ENDIF + + ELSE IF (orientation .EQ. pfeatureType_TrajectoryProfile ) THEN + + IF (along_dim .EQ. x_dim) THEN + plot_orient = pFeatureType_Trajectory + dsg_as_traj = .TRUE. + ELSE IF (along_dim .EQ. z_dim) THEN + GOTO 5000 + ELSE + message = 'PLOT/ALONG='//ww_dim_name(along_dim)// + . ' : TrajectoryProfile data'// + . ' can only be plotted along XY or Z' + ENDIF + + ELSE IF (orientation .EQ. pfeatureType_TimeseriesProfile ) THEN + + IF (along_dim .EQ. x_dim) THEN + plot_orient = pFeatureType_Point + dsg_as_traj = .TRUE. + ELSE IF (along_dim .EQ. z_dim) THEN + GOTO 5000 + ELSE IF (along_dim .EQ. t_dim) THEN + plot_orient = pfeatureType_Timeseries + dsg_as_time = .TRUE. + ELSE + message = 'PLOT/ALONG='//ww_dim_name(along_dim)// + . ' : pfeatureType_TimeseriesProfile data'// + . ' can only be plotted along XY, Z, or T' + ENDIF + + ENDIF + + 4000 CONTINUE + + its_traj = plot_orient.EQ.pfeatureType_Trajectory .OR. + . plot_orient.EQ.pfeatureType_Point .OR. + . dsg_as_traj + + IF (message .NE. 'none') THEN + CALL ERRMSG (ferr_invalid_command, status, message, *5000) + ENDIF + + 5000 RETURN + END diff --git a/fer/plt/dsg_traj_alongxy_setup.F b/fer/plt/dsg_traj_alongxy_setup.F new file mode 100644 index 000000000..fc77b1dac --- /dev/null +++ b/fer/plt/dsg_traj_alongxy_setup.F @@ -0,0 +1,923 @@ + SUBROUTINE DSG_TRAJ_ALONGXY_SETUP ( dset, + . overlay, no_labels, symbol, sym_size, color, + . color1, use_line, do_dash, dashstyle, nokey, + . only_val, skipsym, mv, cx, indep_dat, dep_dat, status ) + +* +* +* This software was developed by the Thermal Modeling and Analysis +* Project(TMAP) of the National Oceanographic and Atmospheric +* Administration's (NOAA) Pacific Marine Environmental Lab(PMEL), +* hereafter referred to as NOAA/PMEL/TMAP. +* +* Access and use of this software shall impose the following +* obligations and understandings on the user. The user is granted the +* right, without any fee or cost, to use, copy, modify, alter, enhance +* and distribute this software, and any derivative works thereof, and +* its supporting documentation for any purpose whatsoever, provided +* that this entire notice appears in all copies of the software, +* derivative works and supporting documentation. Further, the user +* agrees to credit NOAA/PMEL/TMAP in any publications that result from +* the use of this software or in any product that includes this +* software. The names TMAP, NOAA and/or PMEL, however, may not be used +* in any advertising or publicity to endorse or promote any products +* or commercial entity unless specific written permission is obtained +* from NOAA/PMEL/TMAP. The user also understands that NOAA/PMEL/TMAP +* is not obligated to provide the user with any support, consulting, +* training or assistance of any kind with regard to the use, operation +* and performance of this software nor to provide the user with any +* updates, revisions, new versions or "bug fixes". +* +* THIS SOFTWARE IS PROVIDED BY NOAA/PMEL/TMAP "AS IS" AND ANY EXPRESS +* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +* ARE DISCLAIMED. IN NO EVENT SHALL NOAA/PMEL/TMAP BE LIABLE FOR ANY SPECIAL, +* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER +* RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF +* CONTRACT, NEGLIGENCE OR OTHER TORTUOUS ACTION, ARISING OUT OF OR IN +* CONNECTION WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE. +* +* +* set up for line plot by loading data and title and defining axis formats + +* programmer - Ansley Manke +* NOAA/PMEL, Seattle, WA - SDIG +* +* based on dsg_traj_plot_set_up, Set up to draw trajectory plots from a dsg trajectoryProfile dataset +* V762 *acm* 9/20 +* v763 *acm* 9/20 Plots of id-variables in trajectory-profile, timeseries-profile data +* V763 *acm* 9/20 Also plot/along=xy for timeseries or profile data to make a map. + + include 'tmap_dims.parm' + include 'xtm_grid.cmn_text' + external xgt_grid_data + include 'xunits.cmn_text' + include 'ferret.parm' + include 'plot_setup.parm' + include 'errmsg.parm' + include 'slash.parm' + include 'xplot_setup.cmn' + include 'xcontext.cmn' + include 'xdsg_context.cmn' + include 'xvariables.cmn' + include 'xdset_info.cmn_text' + include 'xdyn_linemem.cmn_text' + include 'xprog_state.cmn' + include 'xplot_state.cmn' + include 'xtext_info.cmn' + include 'axis_inc.decl' ! with axis lengths + include 'AXIS.INC' ! with axis lengths + include 'taxis_inc.decl' + include 'TAXIS.INC' + include 'PPLDAT.INC' + include 'switch_inc.decl' + include 'SWITCH.INC' +#include "tmap_dset.parm" + +* calling argument declarations: + LOGICAL overlay, no_labels, use_line, + . nokey + INTEGER dset, symbol, color, color1, do_dash, skipsym, + . mv, cx, status + REAL sym_size, dashstyle(*), only_val + REAL indep_dat(*), dep_dat(*) ! dynamic space allocation + +* V500 *kob* 3/99- up VAR_CODE and VAR_UNITS to 64 chars +* internal variable declarations: + + LOGICAL NO_LINE_RANGE, TM_ITSA_DSG_RAGGED, + . versus, all_1_ind, all_1_dep, formatted, use_nice, + . this_no_range(max_line_on_plot), indep_is_log, set_lev, + . its_dsg, its_cmpnd, mask_applied, user_hlim, is_dsg, + . is_cmpnd, dirs(4), tsprof, its_traj, dsg_as_time, + . dsg_as_traj, as_traj, color_by_station + . + CHARACTER LEFINT*8, LEFT_REAL*16, FULL_VAR_TITLE*2040, + . indep_ax*1, dep_ax*1, tstyle*3, + . buff1*40, buff2*16, t1_date*14, tref*14, fmt*4 + INTEGER TM_LENSTR1, CGRID_SIZE, CGRID_AXIS, DSG_WHATS_IT, + . TM_DSG_DSET_FROM_GRID, TM_DSG_NFEATURES, TM_DSG_OBS_GRID, + . slen, ito, indep_lab, dep_lab, nline_in_mem, + . indep_dim, npts, dep_dim, npts2, ndv, var1, + . grid, ipl, ndim, dims(6), + . slen2, slen3, tax, plot_mem_used, i, + . pxlim, pylim, phlim, pvlim, pindeplim, pdeplim, + . dep_axtyp, indep_axtyp, the_taxis, ribbon_var, + . gap_var, iiunits, iaxis, hblk1, mvh_temp, + . nkey_entries, itmp, igrd_save, cal_id_1, + . icode, dset_dsg, nfeatures, obsdimlen, fline, + . idim, coord_lm, ivar, dim(nferdims), numcx, + . dep_nok, indep_nok, line_symbol, nmasked, nftrsets, + . row_size_lm, orientation, ifeature, imask, iftr, + . lm_index, base, iline, ic, irow_f, irow_l, nobs, + . irow_start, irow_end + + + INTEGER iunits, junits + CHARACTER*1 axdir(nferdims), ax1 + + REAL val1, dt_min, lo, hi, dep_len, + . ind_len, first, last, hlen, dep_lo, + . dep_hi, indep_lo, indep_hi, pad + REAL ind_min, ind_max, dep_min, dep_max, + . delta, tmp, labsize, bad, psum, pcnt + REAL*4 rbad + +* local parameter declarations: + LOGICAL range_rqd + PARAMETER ( range_rqd = .TRUE. ) + DATA axdir / 'X', 'Y', 'Z', 'T', 'E', 'F' / + + INTEGER, DIMENSION(:), ALLOCATABLE :: cmpnd_index + REAL, DIMENSION(:), ALLOCATABLE :: coord + LOGICAL, DIMENSION(:), ALLOCATABLE :: process_feature + LOGICAL, DIMENSION(:), ALLOCATABLE :: process_obs + +* signal that plot set-up has begun + IF ( mode_diagnostic ) CALL SPLIT_LIST(pttmode_ops, ttout_lun, + . 'setting up plot', 15) + +* initialize + + ind_max = arbitrary_small_val8 + ind_min = arbitrary_large_val8 + versus = .TRUE. + indep_dim = plot_axis(1) + dep_dim = plot_axis(2) + nline_in_mem = 0 + phlim = qual_given( slash_hlimits ) + pvlim = qual_given( slash_vlimits ) + pxlim = qual_given( slash_xlimits ) ! deprecated + pylim = qual_given( slash_ylimits ) ! deprecated + dep_axtyp = 1 + indep_axtyp = 1 + indep_is_log = .FALSE. + mod_vs_x = .FALSE. + mod_vs_y = .FALSE. + nkey_entries = 1 + itmp = mnormal + icode = 0 + ribbon_var = 1 + hlen = xlen + line_symbol = symbol + + ! dummy values + cal_id_1 = 0 + t1_date = ' ' + tref = ' ' + +* For DSG dataset get number of features and obsdimlen + + its_dsg = .TRUE. + grid = cx_grid(cx) + color_by_station = DSG_WHATS_IT(grid) .EQ. pdsg_fs_dim + + grid = dsg_xlate_grid(dset) + + CALL TM_DSG_FACTS( grid, idim, obsdimlen, fline, its_dsg, its_cmpnd ) + + IF (its_cmpnd) THEN + CALL TM_DSG_FTRSET_FACTS( grid, dset, orientation, + . nfeatures, nftrsets, is_dsg, is_cmpnd, status ) + ELSE + nftrsets = 0 + nfeatures = TM_DSG_NFEATURES( grid ) + orientation = dsg_orientation(dset) + ENDIF + + tsprof = (orientation .EQ. pfeatureType_TimeseriesProfile) + + CALL DSG_PLTALONG_SETUP (dset, idim, its_traj, dsg_as_traj, + . dsg_as_time, status) + as_traj = orientation.EQ.pfeatureType_TrajectoryProfile .AND. idim.EQ.x_dim + + + ul_dolab(z_dim) = .TRUE. + + ALLOCATE (process_obs(obsdimlen )) + ALLOCATE (process_feature(nfeatures )) + ALLOCATE (cmpnd_index(nfeatures)) + ALLOCATE (coord(nfeatures)) + +* check for improper data supplied + IF ( overlay ) THEN + CALL GET_CX_DIMS( cx, ndim, dims ) + IF ( twodee_on ) indep_dim = dims(1) + IF ( ndim.NE.1 .OR. dims(1).NE.indep_dim ) GOTO 5120 + ELSE + ndim = nplot_axis + IF ( ndim .GT. 1 ) GOTO 5110 + ENDIF + +* Overlay plots, keep the setting of IAUTOT that was set on the underlay plot. + + IF (overlay) THEN + iautot_save = iautot + iautot = 0 + ENDIF + + dep_ax = 'Y' + indep_ax = 'X' + dep_lab = ppl_ylab + indep_lab = ppl_xlab + ind_len = xlen + dep_len = ylen + IF ( pxlim.GT.0 ) THEN + pindeplim = pxlim + ELSE + pindeplim = phlim + ENDIF + IF ( pylim.GT.0 ) THEN + pdeplim = pylim + ELSE + pdeplim = pvlim + ENDIF + +* do not re-draw the main label, logo, other labels from the underlay plot + IF (overlay) THEN + nlabs_on = 0 + CALL PPLCMD ( from, line, 0, 'LABS ', 1, 1 ) + ENDIF + + IF ( no_labels ) CALL PPLCMD ( from, line, 0, 'LABS ', 1, 1 ) + + var1 = 1 ! The variable sent in is the color-by variable. + ndv = 1 ! # of dependent variables + + +* This will cause the x-y context limits to be set up + DO idim = 1, 4 + dirs(idim) = .TRUE. + ENDDO + dirs(z_dim) = .FALSE. + CALL MAKE_DSG_FEATURE_MASK_DIRS(dset, cx, process_feature, nfeatures, dirs) + + IF (is_cmpnd) THEN + lm_index = dsg_loaded_lm(dsg_index_var(dset)) + + DO i = 1, nfeatures + cmpnd_index(i) = dsg_linemem(lm_index)%ptr(i) + 1 + ENDDO + ELSE + DO i = 1, nfeatures + cmpnd_index(i) = i + ENDDO + ENDIF + + row_size_lm = dsg_loaded_lm(dsg_row_size_var(dset)) + + +* * * * * * * - - - INDEPENDENT AXIS DATA - - - * * * * * * * * * + +* Get the LONGITUDE coordinate data + + idim = 1 + iaxis = grid_line(idim, dsg_xlate_grid(dset)) + iunits = line_unit_code ( iaxis ) + IF (iunits .NE. 4) GOTO 5130 ! must be longitude coordinates + + ivar = dsg_coord_var(idim,dset) + coord_lm = dsg_loaded_lm(ivar) ! line memory table indices + + DO iftr = 1, nfeatures + IF (tsprof) THEN + coord(iftr) = dsg_linemem(coord_lm)%ptr(cmpnd_index(iftr)) + ELSE + coord(iftr) = dsg_linemem(coord_lm)%ptr(iftr) + ENDIF + ENDDO + + npts = nfeatures + +* Check for longitudes crossing the dateline. Looks only at location data in features +* included in any feature-level mask, which is applied in PLOT_DSG_APPLY_OBS_MASK + + user_hlim = pindeplim .GT. 0 +c CALL PLOT_DSG_CHECK_LON360 (dset, cx, indep_dat, bad_val4, +c . npts, nfeatures, ind_min, ind_max, user_hlim) + +* Note process_fteset mask masks out only the trajectories masked out by a traj-mask +* not the region spec + + nmasked = nfeatures + imask = 0 + + DO i = 1, nfeatures + IF (process_feature(i)) THEN + imask = imask+1 + indep_dat(imask) = coord(i) + ENDIF + ENDDO + + nmasked = imask + + IF (nmasked .EQ. 0) THEN + cxdsg_empty_set = .TRUE. + + indep_nok = 0 + nmasked = 2 + sym_size = 0.001 + nokey = .TRUE. + qual_given(slash_plot_nokey) = 1 + nokey = .TRUE. + + delta = unspecified_val8 + ind_min = cx_lo_ww(x_dim, cx) + ind_max = cx_hi_ww(x_dim, cx) + + ! Check against dsg lon coord range usin MODSCAT as in MAKE_DSG_OBS_MASK + + indep_dat(1) = bad_val4 + indep_dat(2) = bad_val4 + ELSE + mask_applied = nmasked .LT. nfeatures + CALL MINMAX( indep_dat, nmasked, bad_val4, ind_min, ind_max, indep_nok ) + +* ... force axis scaling if the data has no range + all_1_ind = NO_LINE_RANGE( indep_dat, nmasked, mr_bad_data(mv), val1 ) + delta = unspecified_val8 + + IF ( pindeplim .GT. 0 ) THEN + IF ( pxlim.GT.0 .AND. .NOT.denig_xylim_msg_done ) THEN + CALL WARN( '/XLIMITS and /YLIMITS are deprecated.') + CALL WARN( 'Use /HLIMITS and /VLIMITS instead.') + denig_xylim_msg_done = .TRUE. + ENDIF + + CALL EQUAL_RANGE( + . cmnd_buff(qual_start(pindeplim):qual_end(pindeplim)), + . indep_dim, ind_min, ind_max, delta, + . formatted, range_rqd, cal_id_1, status ) + + IF ( status .NE. ferr_OK ) GOTO 5000 + + ELSEIF ( all_1_ind .AND. .NOT.overlay ) THEN + IF ( val1 .EQ. mr_bad_data(mv) ) THEN + delta = 1. + ind_min = 0.D0 + ind_max = 1.D0 + ELSE + delta = 1. + ind_min = val1 - delta + ind_max = val1 + delta + ENDIF + ELSE + CALL MINMAX( indep_dat, nmasked, bad_val4, lo, hi, indep_nok ) + IF (indep_nok .EQ. 0) THEN ! all missing + ind_min = -1. + ind_max = 1. + ELSE + ind_min = lo + ind_max = hi + ENDIF + ENDIF + + ENDIF ! cxdsg_empty_set + +* Set up independent axis for plotting + + IF (.NOT.overlay) THEN + + IF (.NOT. all_1_ind) THEN + pad = (ind_max - ind_min)* 0.01 + ind_min = ind_min - pad + ind_max = ind_max + pad + ENDIF + + CALL AXIS_ENDS(indep_ax,indep_dim,grid,ind_min,ind_max, + . delta, indep_is_log, indep_axtyp, versus, status) +* This routine checks the units and the setting for formatted lon/lat axes +* If the formatting has been turned off, resets iunits and flag mod_vs_x. + + CALL GET_AXIS_FORMAT( ind_min, ind_max, delta, fmt, use_nice ) + IF (cxdsg_empty_set) use_nice = .FALSE. + IF (use_nice) THEN + ppl_buff = 'XFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LONE'''')' + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + CALL LON_LAT_FMT (x_dim, 'X') + mod_vs_x = .TRUE. + ENDIF + + CALL PPLCMD ( from, line, 0, 'YFOR 0', 1, 1 ) ! for now, anyway + + ELSE ! overlay + CALL PPLCMD ( from, line, 0, indep_ax//'FOR 0', 1, 1 ) + + mod_vs_x = .TRUE. + mod_vs_y = .FALSE. + + ENDIF + +* set up PLOT5 to ignore bad data flag + ppl_buff = ' ' + WRITE ( ppl_buff, 3005 ) bad_val4, indep_ax + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1) + + 3005 FORMAT ( 'LIMITS ',G15.8,1X,A1,'EQ' ) + + ind_min = val1 ! used in the case of only 1 valid data + +* * * * * * * - - - LATITUDE DATA - - - * * * * * * * * * + + IF (cxdsg_empty_set) THEN + npts2 = 2 + dep_nok = 0 + dep_min = cx_lo_ww(y_dim, cx) + dep_max = cx_hi_ww(y_dim, cx) + dep_dat(1) = bad_val4 + dep_dat(2) = bad_val4 + + ELSE + + idim = 2 + iaxis = grid_line(idim, dsg_xlate_grid(dset)) + iunits = line_unit_code ( iaxis ) + IF (iunits .NE. 4) GOTO 5130 ! must be latitude coordinates + + ivar = dsg_coord_var(idim,dset) + coord_lm = dsg_loaded_lm(ivar) ! line memory table indices + + DO iftr = 1, nfeatures + IF (tsprof) THEN + coord(iftr) = dsg_linemem(coord_lm)%ptr(cmpnd_index(iftr)) + ELSE + coord(iftr) = dsg_linemem(coord_lm)%ptr(iftr) + ENDIF + ENDDO + + all_1_dep = .TRUE. + only_val = bad_val4 + + imask = 0 + DO i = 1, nfeatures + IF (.NOT.process_feature(i)) CYCLE + imask = imask + 1 + dep_dat(imask) = coord(i) + ENDDO + +* Set up the dependent axis. + + all_1_dep = NO_LINE_RANGE( dep_dat, nmasked, bad_val4, val1 ) + CALL MINMAX( dep_dat, nmasked, bad_val4, dep_lo, dep_hi, dep_nok ) + IF (dep_nok .EQ. 0) THEN ! all missing + dep_min = -1 + dep_max = 1 + ELSE + dep_min = dep_lo + dep_max = dep_hi + IF (all_1_dep) THEN + delta = 1. + dep_min = dep_lo - delta + dep_max = dep_hi + delta + ENDIF + ENDIF + IF (.NOT.all_1_dep .AND. dep_nok.NE.0) THEN + pad = (dep_max - dep_min)* 0.01 + dep_min = dep_min - pad + dep_max = dep_max + pad + ENDIF + + ENDIF ! cxdsg_empty_set + + + IF (.NOT.overlay) THEN + + delta = unspecified_val8 + CALL AXIS_ENDS(dep_ax, dep_dim, grid, dep_min, dep_max, + . delta, indep_is_log, dep_axtyp, versus, status) + + CALL GET_AXIS_FORMAT( dep_min, dep_max, delta, + . fmt, use_nice ) + IF (cxdsg_empty_set) use_nice = .FALSE. + IF (use_nice) THEN + ppl_buff = 'YFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LAT'''')' + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + ENDIF + ENDIF ! not overlay + + IF ( all_1_dep ) only_val = val1 + +* pass the data to PLOT+ + + WRITE ( ppl_buff, 3005 ) bad_val4, dep_ax + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + + IF (overlay) SOVER = .true. + icode = 0 + npts2 = nmasked + IF (indep_nok.EQ.1 .OR. dep_nok.EQ.1) THEN + dep_dat(1) = dep_lo + indep_dat(1) = ind_min + npts2 = 1 + indep_nok = 1 + IF (line_symbol.EQ.qual_off .OR. + . line_symbol.EQ.unspecified_int4) + . line_symbol = qual_on ! draw point with a symbol + ENDIF + + +* set PLOT+ memory required, longitudes, latitudes, color var, gaps var + + plot_mem_used = 2* 4* npts2 + + CALL PPLLDX_envelope(icode,indep_dat,dep_dat,npts2, + . t1_date, tref, dt_min, plot_mem_used) + +* ... increment number of lines on plot + nline_on = nline_on + 1 + nline_in_mem = nline_in_mem + 1 + CALL LINE_STYLE(line_symbol, sym_size, skipsym, color, color1, use_line, + . do_dash, dashstyle, nline_in_mem, nline_on) + +* * * * * * * - - - LOAD DATA - - - * * * * * * * * * + +* load color-by data into dep_dat + + grid = cx_grid( cx ) + +* set up a color-by variable + + IF (cxdsg_empty_set) THEN +* ... The dep_dat setup from above remains in place. + need_histo = .FALSE. + dep_min = -1. + dep_max = 1. + dep_dat(1) = 1. + dep_dat(2) = 1. + + CALL RIBBON_PLOTKEY_SETUP (overlay, status) + ELSE + + +* ... For a dsg dataset, they can ask to color by an observed variable, or a traj-level +* variablesuch as an ID. That variable can be a string. Either way use it to +* create colorbar label strings and then set values to integer trajectory- +* number for drawing the trajectories. + +* determine 6D shape of context, get size of color var + CALL GET_CX_DIMS_ZERO( cx, ndim, dim ) + dep_dim = dim(1) + npts2 = CGRID_SIZE(cx) + + CALL EXTRACT_LINE ( cx, + . memry(mv)%ptr, + . mv, + . dep_dat, + . dep_dim, + . ndim, + . npts2, + . status ) + + IF (color_by_station) THEN + +* Set up to label the plot key + CALL SET_DSG_ID_LEV + . (dset, cx, nftrsets, dep_dat, mr_type(mv), changed_key) + set_lev = .TRUE. + +* Set the cx_plot limits from the command context so the plot labels will be correct + CALL DSG_ID_CX (cx, nftrsets) + +* Put the trajectory-number into dep_dat + + DO i = 1, nfeatures + dep_dat(i) = cmpnd_index(i) + ENDDO + + npts2 = nfeatures + + ENDIF + + IF ( status .NE. ferr_ok ) GOTO 5000 + +* Issue a note about averaging to color-by one value per profile + + IF (dsg_as_traj .AND. npts2.EQ.obsdimlen) THEN + + IF (cx_lo_ww(t_dim, cx_cmnd).EQ.unspecified_val8 .OR. + . cx_lo_ww(t_dim, cx_cmnd).NE.cx_hi_ww(t_dim, cx_cmnd)) THEN + + + IF (orientation.EQ.t_dim) THEN + CALL WARN ( + . 'PLOT/ALONG= with /T=LO:HI colors the plot with AVE of timeseries data in that range') + ELSE + CALL WARN ( + . 'PLOT/ALONG= with /Z=LO:HI colors the plot with AVE of profile data in that range') + ENDIF + + ENDIF + + ENDIF + +* if color-by trajectory-id variable or timeseries-station variable, +* use the index variable + + +* default color key + CALL RIBBON_PLOTKEY_SETUP (overlay, status) + IF (status .NE. ferr_ok) GOTO 5000 + + set_lev = .FALSE. ! will auto-set levels + bad = mr_bad_data(mv) + +* ... Now get the data values to use in the color-by-variable plot along +* the trajectories. + + IF (npts2 .EQ. nfeatures) THEN + + + imask = 0 + DO ifeature = 1, nfeatures + + IF (.NOT.process_feature(ifeature)) CYCLE + + imask = imask + 1 + dep_dat(imask) = dep_dat(ifeature) + + ENDDO + + ENDIF + + IF (npts2 .EQ. obsdimlen) THEN + + iline = 0 + irow_f = 1 + irow_l = 0 + imask = 0 + + DO ifeature = 1, nfeatures + nobs = dsg_linemem(row_size_lm)%ptr(ifeature) ! feature length + IF (nobs .EQ. 0) GOTO 222 + irow_l = irow_f + nobs - 1 + IF (.NOT.process_feature(ifeature)) GOTO 222 + iline = iline + 1 + +* ... get observation-level mask for this feature + + base = irow_f - 1 + CALL MAKE_DSG_OBS_MASK(dset, cx, ifeature, base, + . process_obs(irow_f), nobs) + +* Find first and last points for this feature applying the user's context +* This just masks the range of the dependent variable. Still need to apply the +* process_obs mask to the dependent variable when plotting. + + irow_start = 0 + irow_end = 0 + + DO ic = irow_f, irow_l + IF (.NOT.process_obs(ic) .AND. irow_start.EQ. 0) THEN + CYCLE + ELSE + irow_start = ic + EXIT ! from loop + ENDIF + ENDDO + + DO ic = irow_l, irow_f, -1 + IF (.NOT.process_obs(ic) .AND. irow_end.EQ. 0) THEN + CYCLE + ELSE + irow_end = ic + EXIT ! from loop + ENDIF + ENDDO + + imask = imask + 1 + IF (irow_start.EQ.0 .OR. irow_end.EQ.0) THEN + dep_dat(imask) = bad_val4 + ELSE + + psum = 0. + pcnt = 0. + DO i = irow_start, irow_end + + IF (dep_dat(i) .NE. bad) THEN + psum = psum + dep_dat(i) + pcnt = pcnt + 1. + ENDIF + + ENDDO + + dep_dat(imask) = bad_val4 + IF (psum .GT. 0.) dep_dat(imask) = psum / pcnt + ENDIF + +* ... prepare for the next feature + + 222 irow_f = irow_l + 1 + + ENDDO + + ENDIF ! npts2 .EQ. obsdimlen + + + +* Compute the mean and standard dev for ribbon-color variable, + CALL MINMAX( dep_dat, nmasked, bad_val4, dep_min, dep_max, dep_nok ) + +* Compute the mean and standard dev for ribbon-color variable, +* needed for computing color levels. Results stored in PPLUS common. +* If histogram-based levels are requested, compute the +* histogram bins. +* send REAL*4 rbad to compare with lev_max, lev_min inside compute_mnstd +* If the user set some /LEVELS then compute_mnstd just returns + + IF (need_histo) THEN +* create temporary buffer to contain workspace + CALL CREATE_TEMP_MEM_VAR( cx, mvh_temp, status ) + IF ( status .NE. ferr_ok ) RETURN + plot_mem_used = plot_mem_used + 2* nmasked + CALL COMPUTE_HISTO_BINS (dep_dat, memry(mvh_temp)%ptr, + . mr_bad_data(mv), nmasked, status) + +* ... clean up temporary variable + CALL DELETE_VARIABLE( mvh_temp ) + ELSE IF (.NOT.set_lev) THEN + rbad = bad_val4 + CALL COMPUTE_MNSTD (dep_dat, mr_bad_data(mv), need_std, nmasked, rbad, status) + + IF (need_std .AND. status.NE.ferr_ok) THEN +* ... set up for automatic levels + CALL PPLCMD ( from, line, 0, 'LEV,()', 1, 1 ) + CALL USE_LINEAR_LEVELS + status = ferr_ok + ENDIF + + ENDIF + + ENDIF ! cxdsg_empty_set +c---------------------- + +* pass the data to PLOT+ + + WRITE ( ppl_buff, 3005 ) mr_bad_data( mv ), dep_ax + + + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + + IF (overlay) SOVER = .true. + +* See pplldx comments. When loading the color-by variable lon vertical time axis +* this forces the right combination of settings. + + icode = 2 + IF (indep_nok .EQ. 1) THEN + dep_dat(1) = dep_min + npts = 1 + ENDIF + CALL PPLLDX_envelope(icode,indep_dat,dep_dat,nmasked, + . t1_date, tref, dt_min, plot_mem_used) + +* ... increment number of lines on plot + nline_on = nline_on + 1 + nline_in_mem = nline_in_mem + 1 + +* Create a gaps-variable. Any masking is applied as the gaps are marked. + + IF (as_traj) THEN + +* first un-masked trajectory, is trajectory # ic + + irow_end = 0 + DO i = 1, nfeatures + IF (process_feature(i) ) THEN + ic = cmpnd_index(i) + EXIT + ENDIF + ENDDO + + imask = 0 + DO i = 1, nfeatures + imask = imask + 1 + dep_dat(imask) = 0. + IF (cmpnd_index(i).NE.ic .AND. imask.GT.1) THEN + dep_dat(imask-1) = 1. + ic = cmpnd_index(i) + ENDIF + ENDDO + + + CALL PPLLDX_envelope(icode, indep_dat, dep_dat, nmasked, + . t1_date, tref, dt_min, plot_mem_used) + + nline_on = nline_on + 1 + nline_in_mem = nline_in_mem + 1 + + ppl_buff = ' ' + gap_var = 4 + WRITE ( ppl_buff, 3007 ) gap_var + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1) + 3007 FORMAT ( 'GAPLOC ',I5 ) + + ENDIF + +* axis scaling and formatting apply possible user settings +* ... independent axis + + IF ( .NOT.overlay ) THEN + +* ... dependent axis scaling +* force axis scaling if the data has no range + + IF ( pdeplim .GT. 0 ) THEN + IF (pylim.GT.0 .AND. .NOT.denig_xylim_msg_done) THEN + CALL WARN( '/XLIMITS and /YLIMITS are deprecated.') + CALL WARN( 'Use /HLIMITS and /VLIMITS instead.') + denig_xylim_msg_done = .TRUE. + ENDIF + + CALL EQUAL_RANGE( + . cmnd_buff(qual_start(pdeplim):qual_end(pdeplim)), + . dep_dim, dep_min, dep_max, delta, + . formatted, range_rqd, cal_id_1, status ) + IF ( status .NE. ferr_OK ) GOTO 5000 + + CALL AXIS_ENDS( dep_ax, the_taxis, grid, dep_min, dep_max, + . delta, indep_is_log, dep_axtyp, versus, status ) + + ELSEIF ( all_1_dep ) THEN + IF (only_val .EQ. bad_val4) val1 = 0.0 ! 10/99 + delta = 1. + IF (val1 .NE. 0) delta = 0.1* val1 + + CALL AXIS_ENDS( dep_ax, the_taxis, grid,val1-delta,val1+delta,delta, + . indep_is_log, dep_axtyp, versus, status ) + + ENDIF + + ax1 = axdir(indep_dim) + CALL PPLCMD ( from, line, 0, 'SET AX_HORIZ '//ax1, 1, 1 ) + ENDIF + +* TITLES +* Main plot title. The main plot label is the the ribbon color variable +* + hlen = xlen ! single-precision xlen from PPLUS common -> double prec. var + + IF ( .NOT.no_labels ) + . CALL LINE_PLOT_LABELS (1, nkey_entries, ndv, cx, + . this_no_range, overlay, versus, nokey, .FALSE., + . tstyle, cal_id_1, ribbon_var, indep_lab, dep_lab, + . ind_min, dep_len, hlen, nfeatures, dset, + . dsg_feature_var(dset) ) + +* set flag indicating a 2D plot is on the screen, and now that things are set up, +* tell any subsequent plots that it's an XY plot + IF (.NOT.overlay) twodee_on = .TRUE. + + plot_axis(1) = 1 + plot_axis(2) = 2 + +* successful completion + + status = ferr_ok + IF (overlay) iautot = iautot_save + + +* not trajectory after all, will just do a ribbon-plot. + + 4000 CONTINUE + + DEALLOCATE (process_obs) + DEALLOCATE (cmpnd_index) + DEALLOCATE (process_feature) + DEALLOCATE (coord) + +* no need to unwind things at the end of xeq_plot. + IF (color_by_station) grid_is_dsg = .FALSE. + + + RETURN + +* error exit + 5000 CALL PPLCMD ( from, line, 0, 'NLINES', 1, 1 ) ! wipe buffers clean + IF (itmp .NE. mnormal) CALL TM_DEALLO_DYN_LINE( itmp ) + GOTO 4000 + + 5110 dep_ax = LEFINT( ndim, slen ) + CALL ERRMSG( ferr_dim_underspec, status, + . 'specified data is not a line'//pCR// + . ' - its a '//dep_ax(:slen)//'D region: "'// + . cmnd_buff(:len_cmnd)//'"', *5000 ) + 5120 CALL ERRMSG( ferr_dim_underspec, status, + . 'overlay is on a different axis'//pCR// + . '"'//cmnd_buff(:len_cmnd)//'"', *5000 ) + + 5130 CALL ERRMSG( ferr_dim_underspec, status, + . 'longitude,latitude coordinate variables not found'//pCR// + . 'Is this a trajectory dataset?'//pCR// + . '"'//cmnd_buff(:len_cmnd)//'"', *5000 ) + + 5150 CALL ERRMSG( ferr_data_type, status, + . 'Cannot plot data type of string on the observation axis'//pCR// + . '"'//cmnd_buff(:len_cmnd)//'"', *5000 ) + + 5200 CALL ERRMSG( ferr_grid_definition, status, + . 'Data grid is not a DSG grid ', + . *5000 ) + + END diff --git a/fer/plt/plot_dsg_data_set_up.F b/fer/plt/plot_dsg_data_set_up.F index 6ad2f403d..b463ee26f 100644 --- a/fer/plt/plot_dsg_data_set_up.F +++ b/fer/plt/plot_dsg_data_set_up.F @@ -79,10 +79,10 @@ SUBROUTINE PLOT_DSG_DATA_SET_UP( nfeatures, nobs_total, . obsdimlen, dset, gxlate, igrid, ifeature, coord_lm(4), . ivar, row_size_lm, nobs, irow_f, irow_l, iaxis, ic, . irow_start, irow_end, base, icxvar, lm_index, - . itim_f, ista, ntim, mv, dep_dim, npts2 + . itim_f, ista, ntim, mv, dep_dim, npts2, i INTEGER*8 rqst_n - REAL value, val_last, val_edge, vdel, bad + REAL value, val_last, val_edge, vdel, bad, psum, pcnt CHARACTER buff1, LEFINT*8, buff8*8, buff*256 INTEGER, DIMENSION(:), ALLOCATABLE :: station_index @@ -416,15 +416,8 @@ SUBROUTINE PLOT_DSG_DATA_SET_UP( nfeatures, nobs_total, CALL TM_ALLO_TMP_GRID (igrid, status) IF (status .NE. ferr_ok) GOTO 5000 -c IF (dsg_as_time) THEN ! just get the time axis into new grid -c CALL TM_COPY_GRID (grid, igrid) -c DO idim = 1, nferdims -c IF (idim .NE. the_dim) grid_line(idim, igrid) = 0 -c ENDDO -c ELSE - CALL TM_COPY_GRID_W_LINE_USE (grid, igrid) - CALL TM_DEALLO_DYN_LINE( grid_line(the_dim, grid) ) -c ENDIF + CALL TM_COPY_GRID_W_LINE_USE (grid, igrid) + CALL TM_DEALLO_DYN_LINE( grid_line(the_dim, grid) ) CALL TM_ALLO_TMP_LINE(iaxis, status) CALL GET_LINE_DYNMEM (nobs_total, iaxis, status) @@ -481,10 +474,30 @@ SUBROUTINE PLOT_DSG_DATA_SET_UP( nfeatures, nobs_total, ntim = ntim + 1 value = dsg_linemem(coord_lm(the_dim))%ptr(ifeature) CALL PUT_LINE_COORD (linemem(iaxis)%ptr, ntim, value) - CALL PUT_LINE_COORD (dsg_linemem(cxdsg_lm_tsdat)%ptr, ntim, dep_dat(irow_start) ) + + +* Compute the average data from this feature in the z range given + psum = 0. + pcnt = 0. + DO i = irow_start, irow_end + + IF (dep_dat(i) .NE. bad) THEN + psum = psum + dep_dat(i) + pcnt = pcnt + 1. + ENDIF + + ENDDO + + IF (psum .GT. 0.) THEN + psum = psum / pcnt + ELSE + psum = bad_val4 + ENDIF + + CALL PUT_LINE_COORD (dsg_linemem(cxdsg_lm_tsdat)%ptr, ntim, psum ) IF (irow_start .NE. irow_end .AND. first_warn) THEN CALL WARN ( - . 'PLOT/ALONG= /Z=LO:HI uses the first data value in the Z- range on each profile') + . 'PLOT/ALONG= /Z=LO:HI plots the AVE data in the Z- range on each profile') first_warn = .FALSE. ENDIF ELSE diff --git a/fer/plt/plot_vs_set_up.F b/fer/plt/plot_vs_set_up.F new file mode 100644 index 000000000..8019aac81 --- /dev/null +++ b/fer/plt/plot_vs_set_up.F @@ -0,0 +1,1480 @@ + SUBROUTINE PLOT_VS_SET_UP( + . overlay, transpz, no_labels, + . symbol, sym_size, color, color1, use_line, + . step_inc, do_dash, dashstyle, + . is_logx, is_logy, nokey, addgaps, + . all_1_dep, only_val, skipsym, mv_list, + . cx_list, nmv, indep_dat, dep_dat, status ) + +* +* +* This software was developed by the Thermal Modeling and Analysis +* Project(TMAP) of the National Oceanographic and Atmospheric +* Administration's (NOAA) Pacific Marine Environmental Lab(PMEL), +* hereafter referred to as NOAA/PMEL/TMAP. +* +* Access and use of this software shall impose the following +* obligations and understandings on the user. The user is granted the +* right, without any fee or cost, to use, copy, modify, alter, enhance +* and distribute this software, and any derivative works thereof, and +* its supporting documentation for any purpose whatsoever, provided +* that this entire notice appears in all copies of the software, +* derivative works and supporting documentation. Further, the user +* agrees to credit NOAA/PMEL/TMAP in any publications that result from +* the use of this software or in any product that includes this +* software. The names TMAP, NOAA and/or PMEL, however, may not be used +* in any advertising or publicity to endorse or promote any products +* or commercial entity unless specific written permission is obtained +* from NOAA/PMEL/TMAP. The user also understands that NOAA/PMEL/TMAP +* is not obligated to provide the user with any support, consulting, +* training or assistance of any kind with regard to the use, operation +* and performance of this software nor to provide the user with any +* updates, revisions, new versions or "bug fixes". +* +* THIS SOFTWARE IS PROVIDED BY NOAA/PMEL/TMAP "AS IS" AND ANY EXPRESS +* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +* ARE DISCLAIMED. IN NO EVENT SHALL NOAA/PMEL/TMAP BE LIABLE FOR ANY SPECIAL, +* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER +* RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF +* CONTRACT, NEGLIGENCE OR OTHER TORTUOUS ACTION, ARISING OUT OF OR IN +* CONNECTION WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE. +* +* +* set up for line plot by loading data and title and defining axis formats +* This version handles all of the logic for versus plots, pulling that out of +* plot_set_up.F +* 9/2020 ACM +* v763 *acm* 9/20 Fixing bug with PLOT/VS of single-valued data, negative longitudes. + + IMPLICIT NONE + include 'tmap_dims.parm' + include 'xtm_grid.cmn_text' + external xgt_grid_data + include 'xunits.cmn_text' + include 'ferret.parm' + include 'plot_setup.parm' + include 'errmsg.parm' + include 'slash.parm' + include 'xplot_setup.cmn' + include 'xcontext.cmn' + include 'xvariables.cmn' + include 'xprog_state.cmn' + include 'xplot_state.cmn' + include 'xtext_info.cmn' + include 'axis_inc.decl' ! with axis lengths + include 'AXIS.INC' ! with axis lengths + include 'taxis_inc.decl' + include 'TAXIS.INC' + include 'PPLDAT.INC' + include 'switch_inc.decl' + include 'SWITCH.INC' + include 'xdset_info.cmn_text' + include 'xdsg_context.cmn' + +#include "tmap_dset.parm" + +* calling argument declarations: + LOGICAL overlay, transpz, no_labels, use_line, + . is_logx, is_logy, nokey, all_1_dep, addgaps + INTEGER symbol, color, color1, step_inc, do_dash, skipsym, + . nmv, mv_list( * ), cx_list( * ), status + REAL sym_size, + . dashstyle(*), only_val + REAL indep_dat(*), dep_dat(*) ! dynamic space allocation + +* V500 *kob* 3/99- up VAR_CODE and VAR_UNITS to 64 chars +* internal variable declarations: + + LOGICAL NO_LINE_RANGE, GEOG_LABEL, ITS_FMRC, + . ITSA_TRUEMONTH_AXIS, TM_HAS_STRING, TM_ITSA_DSG, + . versus, flip, time_axis, + . all_1_ind, formatted, + . indep_is_log, dep_is_log, use_nice, + . this_no_range(max_line_on_plot), set_axis, + . timeax_both, do_units, doing_ribbon_color, + . x_is_time, y_is_time, its_traj, its_dsg, its_cmpnd, + . expand_dsg_ok, color_dsg_id, user_hlim, dsg_as_time, + . dsg_as_traj, do_traj, native_traj + + CHARACTER*2040 FULL_VAR_TITLE, KEY_STRING, plot_title + CHARACTER INTERNAL_WHOI_DATE*14, VAR_UNITS*64, + . SECS_TO_DATE*20, LEFINT*8, MERGED_WHOI_DATE*14, + . LEFT_REAL*16, TM_GET_CALENDAR_NAME*32, + . indep_ax*1, dep_ax*1, tstyle*3, xtra_lab*24, + . buff1*40, buff2*16, buff3*8, t1_date*14, tref*14, + . cal_name*32, cal_name_new*32, val_buff*45, fmt*4, + . whoimin*14, whoimax*14, + . whoimin_axlo*14, whoimax_axhi*14 + INTEGER TM_LENSTR1, CX_DIM_LEN, CGRID_SIZE, CGRID_AXIS, + . TM_DSG_DSET_FROM_GRID, TM_DSG_NFEATURES, TM_DSG_NF2FEATURES, + . mv, cx, slen, ito, indep_lab, dep_lab, nline_in_mem, + . indep_dim, npts, dep_dim, npts2, ndv, var1, + . grid, grid1, ipl, ndim, dims(6), + . slen2, slen3, tax, plot_mem_used, nload, i, + . pxlim, pylim, phlim, pvlim,pindeplim, pdeplim, + . dep_axtyp, indep_axtyp, the_taxis, ribbon_var, + . gap_var, iiunits, iaxis, hblk1, mvh_temp, + . nkey_entries, t_axtyp, since_T0, itmp, igrd_save, + . icode, dset, nfeatures, obsdimlen, fline, orientation, + . dtype, fvar, nload_all, numv, nok, line_symbol, + . dsg_fmt_grid, dep_dir, nmasked, dset_dsg, dum1, dum2, + . nftrsets + + INTEGER TM_GET_CALENDAR_ID, cal_id_1, cal_id_2, + . cal_id_old, cal_id_new + INTEGER TM_UNIT_ID, cx_x , cx_y, iunits, junits, axis_x_dir, axis_y_dir + CHARACTER*1 axdir(6), ax1 + + REAL val1, dt_min, lo, hi, dep_len, + . ind_len, first, last, hlen, hrs, + . bc_min, bc_max, step_lo, step_hi, ascale, aoff + REAL TSTEP_TO_SECS, SECS_TO_TSTEP, DATE_TO_TSTEP, + . WHOI2BC, GET_LINE_COORD, ind_min, ind_max, dep_min, dep_max, + . delta, pad, tmp, dum, ww_save, cen, yr + REAL*4 rbad + +* local parameter declarations: + LOGICAL main, norm_labs, pnot_curv, range_rqd + PARAMETER ( main = .TRUE., + . norm_labs = .FALSE., + . pnot_curv = .FALSE., + . range_rqd = .TRUE. ) + DATA axdir / 'X', 'Y', 'Z', 'T', 'E', 'F' / + +* signal that plot set-up has begun + IF ( mode_diagnostic ) CALL SPLIT_LIST(pttmode_ops, ttout_lun, + . 'setting up plot', 15) + + versus = .TRUE. + +* limits error check + IF (nmv .GT. max_line_on_plot) GOTO 5050 + +* initialize + + ind_max = arbitrary_small_val8 + ind_min = arbitrary_large_val8 + cx = cx_list( 1 ) + mv = mv_list( 1 ) + grid1 = cx_grid( cx ) + grid = cx_grid(is_cx(1)) + dset = cx_data_set(cx) + indep_dim = plot_axis(1) + twodee_on = .TRUE. + + numv = nmv ! local copy + line_symbol = symbol + IF (cxdsg_empty_set) THEN + sym_size = 0.001 + line_symbol = 1 ! doesn't matter what + use_line = .FALSE. + ENDIF + +* ragged DSG data? Get number of features and obsdimlen. If it's a +* Trajectory or Point dataset and they did "PLOT var", set that up. + + nfeatures = 0 + its_dsg = .FALSE. + IF (dset .EQ. pdset_irrelevant) THEN + cx = is_cx(nmv) + grid1 = cx_grid( cx ) + dset = cx_data_set(cx) + ENDIF + +* a user-var or expression where the grid is dsg? + + IF (dset .GT. pdset_irrelevant) its_dsg = dsg_ragged(dset) + + dset_dsg = TM_DSG_DSET_FROM_GRID( grid ) + IF (dset_dsg .GT. pdset_irrelevant) THEN + IF (.NOT. its_dsg) dset = dset_dsg + its_dsg = dsg_ragged(dset_dsg) + ENDIF + + IF (its_dsg ) THEN + igrd_save = grid + grid = cx_grid(is_cx(nmv)) + IF (dset_dsg .LT. pdset_irrelevant) grid = dsg_xlate_grid(dset) + + nfeatures = TM_DSG_NFEATURES( grid ) + +* The dataset may be DSG but the variable not. + CALL TM_DSG_FACTS( grid, orientation, obsdimlen, fline, its_dsg, its_cmpnd ) + + nftrsets = TM_DSG_NF2FEATURES (dset) + + native_traj = orientation.EQ.pfeatureType_Trajectory .OR. + . orientation.EQ.pfeatureType_Point + +* Is this a PLOT/along=xy for some other dataset type? +* A bit convoluted: do_traj when plotting DSG as an XY map plot for any reason +* CHECK_PLOT_TRAJ is testing for a scatter plot based on a dsg dataset +* Then the XY map plot is called with a flag its_traj for whether it's +* trajectories, with multiple values, or Points with one value per feature. + + CALL DSG_PLTALONG_SETUP (dset, orientation, do_traj, + . dsg_as_traj, dsg_as_time, status) + IF (status .NE. ferr_ok) GOTO 5000 + IF (dsg_as_time) plot_axis(1) = t_dim + its_traj = do_traj + +* Is this just a regular ribon plot with 3 variables (not a lon/lat traj plot)? +* OR a PLOT/VS var1,var2 from a DSG dataset for var1,var2 not lon/lat? +* This may turn off ribbon_plot + + IF (.NOT.dsg_as_traj .AND. .NOT.dsg_as_time) + . CALL CHECK_PLOT_TRAJ (versus, numv, grid, its_traj) + + IF (.NOT.its_traj) native_traj = .FALSE. + IF (its_traj) ribbon_plot = .TRUE. + +* numv = 1 if we're doing the trajectory ribbon plot automatically with PLOT var +* Still allow doing it the old fashioned way numv = 2 for plot/vs lon, lat +* If it's a trajectory dataset we still might send in say, lon,lat with a map projection +* specified, so allow for that too. + + IF ( (numv.EQ.1 .OR. nmv.EQ.3) .AND. its_traj) THEN + ul_dolab(x_dim) = .FALSE. + ul_dolab(y_dim) = .FALSE. + ul_dolab(t_dim) = .TRUE. + ul_dolab(e_dim) = .TRUE. + ul_nlabs = 2 + + IF (orientation .EQ. pfeatureType_Point) THEN + use_line = .FALSE. + IF (symbol .EQ. unspecified_int4) symbol = 20 + ENDIF + + IF (native_traj) THEN + its_traj = dsg_feature_type(dset) .EQ. pfeatureType_Trajectory ! as opposed to Point + CALL DSG_TRAJ_PLOT_SET_UP ( + . overlay, no_labels, symbol, sym_size, color, color1, + . use_line, do_dash, dashstyle, nokey, its_traj, only_val, + . skipsym, mv_list(1), cx_list(1), numv, indep_dat, dep_dat, status ) + + IF (numv.EQ.3 .AND. status.NE.ferr_ok) GOTO 2000 ! not traj after all, just a ribbon plot + + ELSE IF (do_traj .OR. dsg_as_traj) THEN + CALL DSG_TRAJ_ALONGXY_SETUP( dset, + . overlay, no_labels, symbol, sym_size, color, + . color1, use_line, do_dash, dashstyle, nokey, + . only_val, skipsym, mv_list, cx_list, + . indep_dat, dep_dat, status ) + + ENDIF + + GOTO 1000 + + ENDIF ! its_traj + + grid = igrd_save + + ENDIF ! its_dsg + +* Plot a variable e.g. as a set of timeseries per feature + +* initialize + + 2000 CONTINUE + + dep_max = arbitrary_small_val8 + dep_min = arbitrary_large_val8 + ind_max = arbitrary_small_val8 + ind_min = arbitrary_large_val8 + cx = cx_list( 1 ) + mv = mv_list( 1 ) + grid1 = cx_grid( cx ) + indep_dim = plot_axis(1) + IF (dsg_as_time) THEN + indep_dim = t_dim + time_axis = .TRUE. + cal_id_1 = 1 + ENDIF + nline_in_mem = 0 + phlim = qual_given( slash_hlimits ) + pvlim = qual_given( slash_vlimits ) + pxlim = qual_given( slash_xlimits ) ! deprecated + pylim = qual_given( slash_ylimits ) ! deprecated + cal_id_1 = 0 + cal_id_2 = 0 + cal_name = ' ' + dep_axtyp = 1 + indep_axtyp = 1 + mod_vs_x = .FALSE. + mod_vs_y = .FALSE. + axis_x_dir = no_dim + axis_y_dir = no_dim + adjust_time = .FALSE. + cx_x = cx_list(1) + xtra_lab = ' ' + nkey_entries = numv + itmp = mnormal + do_units = .TRUE. + x_is_time = .FALSE. + y_is_time = .FALSE. + tstyle = ' ' + the_taxis = no_dim + icode = 0 + ribbon_var = 0 + gap_var = 0 + t1_date = ' ' + tref = ' ' + +* set limit on the number of plot keys + IF (.NOT.nokey .AND. numv.GT.max_key_entries) THEN + buff1 = LEFINT( max_key_entries, slen ) + CALL WARN('Too many lines to label. Labeling initial ' + . //buff1(:slen) ) + nkey_entries = max_key_entries + ENDIF + +* determine PLOT+ memory required + plot_mem_used = 0 + DO 10 ipl = 1, numv + 10 plot_mem_used = plot_mem_used + . + 2*step_inc*CGRID_SIZE( cx_list(ipl) ) + +* check for improper data supplied + IF ( num_uvars_in_cmnd .LT. 2 ) GOTO 5100 + +* length of plot array + npts = CGRID_SIZE( cx ) ! (redundantly calc'd below) 5/90 + +* Overlay plots, keep the setting of IAUTOT that was set on the underlay plot. + + IF (overlay) THEN + iautot_save = iautot + iautot = 0 + ENDIF + +* Mark the gap-variable for PLOT/VS + + IF (grid_is_dsg .AND. .NOT.its_traj) addgaps = .TRUE. + + IF (addgaps) THEN + IF (ribbon_plot) THEN + gap_var = 4 + ELSE + IF (num_uvars_in_cmnd .GT. 3) GOTO 5190 + gap_var = 3 + ENDIF + ENDIF + +* set flag indicating a 1D plot is on the screen + IF (.NOT.overlay) onedee_on = .TRUE. + + IF (ribbon_plot) ribbon_var = 3 + +* decide if the picture should be rotated 90 degrees + flip = transpz + + IF ( flip ) THEN + dep_ax = 'X' + indep_ax = 'Y' + dep_lab = ppl_xlab + indep_lab = ppl_ylab + dep_len = xlen + ind_len = ylen + indep_is_log = is_logy + dep_is_log = is_logx + + IF ( pxlim.GT.0 ) THEN + pindeplim = pxlim + ELSE + pindeplim = pvlim + ENDIF + IF ( pylim.GT.0 ) THEN + pdeplim = pylim + ELSE + pdeplim = phlim + ENDIF + ELSE + dep_ax = 'Y' + indep_ax = 'X' + dep_lab = ppl_ylab + indep_lab = ppl_xlab + ind_len = xlen + dep_len = ylen + indep_is_log = is_logx + dep_is_log = is_logy + IF ( pxlim.GT.0 ) THEN + pindeplim = pxlim + ELSE + pindeplim = phlim + ENDIF + IF ( pylim.GT.0 ) THEN + pdeplim = pylim + ELSE + pdeplim = pvlim + ENDIF + ENDIF + +* do not re-draw the main label, logo, other labels from the underlay plot + IF (overlay) THEN + nlabs_on = 0 + CALL PPLCMD ( from, line, 0, 'LABS ', 1, 1 ) + ENDIF + + IF ( no_labels ) CALL PPLCMD ( from, line, 0, 'LABS ', 1, 1 ) + +* * * * * * * - - - INDEPENDENT AXIS DATA - - - * * * * * * * * * + + time_axis = .FALSE. + adjust_time = .TRUE. ! (Says, do not try to further adjust time...) + CALL EXTRACT_LINE ( cx, + . memry(mv)%ptr, + . mv, + . indep_dat, + . indep_dim, + . ndim, + . npts, + . status ) + IF ( status .NE. ferr_ok ) GOTO 5000 + var1 = 2 ! since first var is independent axis + ndv = numv - 1 ! # of dependent variables + use_keys = .FALSE. + plot_title = FULL_VAR_TITLE( cx_list(1), .TRUE., slen ) + +* ... Is this a dsg Profile dataset? Are we doing a vs plot? If so we may need to expand the +* first or second variable to the OBS axis. + + expand_dsg_ok = .TRUE. + IF (its_dsg) THEN ! Get the size of the second variable and decide whether to expand to obs axis. + + npts2 = CGRID_SIZE( cx_list(var1) ) + + expand_dsg_ok = (npts.NE.nftrsets .AND. + . (npts.NE.nfeatures .OR. npts2.NE.nfeatures) ) + dsg_fmt_grid = 0 + + IF (its_cmpnd) orientation = pfeatureType_Profile + IF ( orientation.EQ.pfeatureType_Profile .AND. + . versus .AND. expand_dsg_ok) THEN + + IF (npts .EQ. nfeatures) CALL DSG_OBS_BY_FEATURE_VAR + . (dset, cx, nfeatures, obsdimlen, indep_dat) + npts = obsdimlen + dsg_fmt_grid = dsg_xlate_grid(dset) + ENDIF + IF (dsg_feature_type(dset).EQ.pfeatureType_Point) expand_dsg_ok = .TRUE. + ENDIF + + +* ... Is this a dsg dataset and we're doing a DSG plot (not a native trajectory plot) then +* apply an obs mask. Resets npts. + + IF (grid_is_dsg .AND. .NOT.its_traj .AND.expand_dsg_ok) + . CALL PLOT_DSG_APPLY_OBS_MASK (dset, cx, nfeatures, obsdimlen, + . indep_dat, npts) + +* ... label the independent axis as the first var given + + buff1 = VAR_UNITS(cx_x) + iunits = TM_UNIT_ID(buff1) + +* Does it have time units? If so draw a formatted time axis + x_is_time = .FALSE. + since_T0 = MAX( INDEX(buff1,'since'), INDEX(buff1,'SINCE') ) + slen = TM_LENSTR1(buff1) + IF (since_T0.GT.0 .AND. slen.GT.since_T0+5) THEN + cal_id_1 = 1 + CALL TM_ALLO_DYN_LINE( itmp, status ) + CALL CD_GET_TIME_UNITS ( buff1, cal_id_1, line_units(itmp), + . line_t0(itmp), dum, status ) + IF (status .EQ. ferr_ok) THEN + x_is_time = .TRUE. + the_taxis = indep_dim + iunits = TM_UNIT_ID( line_units(itmp) ) + line_unit_code(itmp) = TM_UNIT_ID( line_units(itmp) ) + line_tunit(itmp) = un_convert(line_unit_code(itmp)) + line_cal_name(itmp) = 'GREGORIAN' + line_direction( itmp ) = 'TI' + + CALL MINMAX( indep_dat, npts, mr_bad_data(mv), lo, hi, nok ) + grid_line(the_taxis, mgrid_buff) = itmp + + igrd_save = cx_grid(cx_list(1)) + cx_grid(cx_list(1)) = mgrid_buff + ww_save = cx_lo_ww(the_taxis,cx_list(1) ) + cx_lo_ww(the_taxis,cx_list(1)) = lo + + ind_min = TSTEP_TO_SECS (mgrid_buff, the_taxis, lo) + ind_max = TSTEP_TO_SECS (mgrid_buff, the_taxis, hi) + + IF (dsg_fmt_grid .GT. 0) THEN + ind_min = TSTEP_TO_SECS (dsg_fmt_grid, t_dim, lo) + ind_max = TSTEP_TO_SECS (dsg_fmt_grid, t_dim, hi) + ENDIF + +* Plotting short time range? If so recompute tref +* axis length in hours from time since BC in seconds + IF (.NOT.overlay) THEN + hrs = ( ind_max - ind_min ) / 3600. + IF (hrs.LT.24) tref = + . MERGED_WHOI_DATE( cx_list(var1), the_taxis, numv, .TRUE. ) + + all_1_ind = ind_min .EQ. ind_max + IF ( all_1_ind ) THEN + ind_min = ind_min - 1. + ind_max = ind_max + 1. + ENDIF + + CALL TAXIS_STYLE( indep_ax, ind_min, ind_max, tstyle, xtra_lab ) +* 11/16 Get the min and max as used by PPLUS. That function uses WHOI formatted time + CALL TPLOT_AXIS_ENDS (ind_min, ind_max, cal_id_1, tstyle) + + IF (.NOT.overlay) CALL AXIS_END_SYMS( indep_ax, + . SECS_TO_TSTEP( mgrid_buff, the_taxis, ind_min ), + . SECS_TO_TSTEP( mgrid_buff, the_taxis, ind_max ) ) + + ENDIF + + ELSE + IF (itmp .NE. mnormal) CALL TM_DEALLO_DYN_LINE( itmp ) + itmp = mnormal + ENDIF + ELSE ! not time axis + + CALL MINMAX( indep_dat, npts, mr_bad_data(mv), lo, hi, nok ) + ind_min = MIN(ind_min, lo) + ind_max = MAX(ind_max, hi) + ENDIF + + slen = TM_LENSTR1(plot_title) + +* Check for lon-360 issues? Crossing the dateline within features, +* when longitude is in the obs direction. + + CALL GEOG_LABEL_VS (buff1, iunits, x_dim, axis_x_dir) + + IF (iunits .EQ. 4) THEN + delta = MAX(1., (ind_max-ind_min)/10.) ! just a rough guess, will be set later + CALL GET_AXIS_FORMAT( ind_min, ind_max, delta, fmt, use_nice ) + + IF (grid_is_dsg .AND. versus .AND. use_nice .AND. + . axis_x_dir.EQ.x_dim .AND. npts.NE.nfeatures) THEN + + user_hlim = pindeplim .GT. 0 + CALL PLOT_DSG_CHECK_LON360 (dset, cx, indep_dat, mr_bad_data(mv), + . npts, nfeatures, dum1, dum2, .TRUE.) + + ENDIF + ENDIF + + IF (.NOT.no_labels .AND. iunits .NE. 4) + . CALL BOX_LABEL( indep_lab, + . plot_title(:slen), + . 0.0, 0.0, 0.6*ind_len, + . dflt_letsize_label*textscale, + . dflt_letsize_label*textscale, + . ppl_centered, + . lab_loc_absolute, lab_loc_absolute ) + + +* ... force axis scaling if the data has no range + all_1_ind = NO_LINE_RANGE( indep_dat,npts,mr_bad_data(mv),val1 ) + delta = unspecified_val8 + + IF ( pindeplim .GT. 0 ) THEN + IF ( pxlim.GT.0 .AND. .NOT.denig_xylim_msg_done ) THEN + CALL WARN( '/XLIMITS and /YLIMITS are deprecated.') + CALL WARN( 'Use /HLIMITS and /VLIMITS instead.') + denig_xylim_msg_done = .TRUE. + ENDIF +* Get calendar name for equal_range (even if not time axis) + tax = grid_line(f_dim,cx_grid(cx_list(var1)) ) + IF (tax .EQ. 0) tax = grid_line(t_dim,cx_grid(cx_list(var1)) ) + cal_name = line_cal_name(tax) + cal_id_1 = TM_GET_CALENDAR_ID ( cal_name ) + + CALL EQUAL_RANGE( + . cmnd_buff(qual_start(pindeplim):qual_end(pindeplim)), + . indep_dim, ind_min, ind_max, delta, + . formatted, range_rqd, cal_id_1, status ) + +* Check for a valid axis if they've asked for a log axis before proceeding + + IF (indep_is_log) THEN + CALL AXIS_ENDS(indep_ax,indep_dim,grid1,ind_min,ind_max, + . delta, indep_is_log, indep_axtyp, versus, status) + IF ( status .NE. ferr_OK ) THEN + first = ind_min + last = ind_max + GOTO 5170 + ENDIF + +* Get requested ends in data units again + CALL EQUAL_RANGE( + . cmnd_buff(qual_start(pindeplim):qual_end(pindeplim)), + . indep_dim, ind_min, ind_max, delta, + . formatted, range_rqd, cal_id_1, status ) + + ENDIF + + IF ( status .NE. ferr_OK ) GOTO 5000 + ELSEIF ( all_1_ind .AND. .NOT.overlay ) THEN + IF ( val1 .EQ. mr_bad_data(mv) ) THEN + delta = 1. + ind_min = 0.D0 + ind_max = 1.D0 + ELSE + delta = 0.1* ABS(val1) + ind_min = val1 - delta + ind_max = val1 + delta + ENDIF + ELSE + CALL MINMAX( indep_dat, npts, mr_bad_data(mv), lo, hi, nok ) + IF (nok .EQ. 0) THEN ! all missing + ind_min = -1. + ind_max = 1. + ElSE + ind_min = lo + ind_max = hi + ENDIF + ENDIF + + IF (.NOT. all_1_ind .AND. .NOT.indep_is_log) THEN + pad = (ind_max - ind_min)* 0.01 + ind_min = ind_min - pad + ind_max = ind_max + pad + ENDIF + +* Set up modulo longitude handling if /vs and variable has units +* of longitude + + cx_x = cx_list(1) + cx_y = cx_list(2) + buff1 = VAR_UNITS(cx_x) + IF (iunits .GE. 0) iunits = TM_UNIT_ID(buff1) + + + IF (.NOT.overlay) THEN + + CALL AXIS_ENDS(indep_ax,indep_dim,grid1,ind_min,ind_max, + . delta, indep_is_log, indep_axtyp, versus, status) + IF ( status .NE. ferr_ok ) THEN + first = ind_min + last = ind_max + GOTO 5170 + ENDIF + +* This routine checks the units and the setting for formatted lon/lat (time?) axes +* If the formatting has been turned off, resets iunits and flag mod_vs_x. +* If the axes are labeled with longitude, latitude, or time, then set ul_dolab to not +* put upper-left labels on the plot with the ranges. + + CALL GEOG_LABEL_VS (buff1, iunits, x_dim, axis_x_dir) + + IF (iunits .EQ. 4) THEN + CALL GET_AXIS_FORMAT( ind_min, ind_max, delta, + . fmt, use_nice ) + IF (use_nice) THEN + plot_axis(1) = 1 + IF (axis_x_dir .EQ. x_dim) THEN + ppl_buff = 'XFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LONE'''')' + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + CALL LON_LAT_FMT (x_dim, 'X') + mod_vs_x = .TRUE. + ul_dolab(x_dim) = .FALSE. + + ELSEIF (axis_x_dir .EQ. y_dim) THEN + ppl_buff = 'XFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LAT'''')' + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + CALL LON_LAT_FMT (y_dim, 'Y') + ul_dolab(y_dim) = .FALSE. + plot_axis(1) = 2 + ELSE + ppl_buff = 'XFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LONE'''')' + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + CALL LON_LAT_FMT (x_dim, 'X') + ul_dolab(x_dim) = .FALSE. + ENDIF + ENDIF + ELSEIF (iunits .LT. 0) THEN ! time axis formatting + CALL PPLCMD ( from, line, 0, 'TIME', 1, 1 ) + CALL PPLCMD ( from, line, 0, 'TAXIS 60,ON', 1, 1 ) + ppl_buff = ' ' + WRITE ( ppl_buff, 3004 ) 'GREGORIAN' + CALL PPLCMD (from, line, 0, ppl_buff, 1, 1) + ul_dolab(t_dim) = .FALSE. + plot_axis(1) = 4 + ELSE + CALL PPLCMD ( from, line, 0, 'XFOR 0', 1, 1 ) + ENDIF + + CALL PPLCMD ( from, line, 0, 'YFOR 0', 1, 1 ) ! for now, anyway + + ELSE ! overlay + CALL PPLCMD ( from, line, 0, indep_ax//'FOR 0', 1, 1 ) + + IF (iunits .EQ. 4) THEN + mod_vs_x = (TM_HAS_STRING(buff1, '_e') .OR. + . TM_HAS_STRING(buff1, 'lon') ) + + ENDIF + + buff1 = VAR_UNITS(cx_y) + junits = TM_UNIT_ID(buff1) + IF (junits .EQ. 4) THEN + mod_vs_y = (TM_HAS_STRING(buff1, '_e') .OR. + . TM_HAS_STRING(buff1, 'lon') ) + ENDIF + + ENDIF + +* set up PLOT5 to ignore bad data flag + ppl_buff = ' ' + WRITE ( ppl_buff, 3005 ) mr_bad_data( mv ), indep_ax + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1) + + + 3004 FORMAT ('CALENDAR ', A7) + 3005 FORMAT ( 'LIMITS ',G15.8,1X,A1,'EQ' ) + +* Send info about gap var for PLOT/VS/GAPLOC + + CALL PPLCMD ( from, line, 0, 'GAPLOC 0', 1, 1) + IF (addgaps) THEN + ppl_buff = ' ' + WRITE ( ppl_buff, 3007 ) gap_var + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1) + ENDIF + 3007 FORMAT ( 'GAPLOC ',I5 ) + +* * * * * * * - - - LOAD DATA - - - * * * * * * * * * +* 11/95: determine T0 reference time for time axis + IF ( time_axis .AND. .NOT.overlay ) THEN + tref = MERGED_WHOI_DATE( cx_list(var1), the_taxis, numv, .FALSE. ) + ELSE + tref = ' ' ! overlays reuse previous tref + ENDIF + +* loop through all the dependent variables loading them into PPLUS + + all_1_dep = .TRUE. + set_axis = .TRUE. + only_val = bad_val4 + nload_all = 0 + + IF (grid_is_dsg .AND. .NOT.its_traj .AND. use_line ) numv = numv + 1 ! fake a gap-var + + DO 200 ipl = var1, numv + IF (grid_is_dsg .AND. .NOT.its_traj .AND. ipl.EQ.gap_var) THEN + +* create the gaps-var. IF npts2 = nfeatures (e.g. a plot of lon vs lat for timeseries stations) +* then the gap-var will be all 0's, e.g. no gaps. + cx = cx_list( 1 ) + IF (npts2 .NE. nfeatures) npts2 = obsdimlen + CALL DSG_OBS_MARK_GAPS (dset, cx, nfeatures, npts2, dep_dat) + status = ferr_ok + + ELSE + mv = mv_list( ipl ) + cx = cx_list( ipl ) + grid = cx_grid( cx ) + dset = cx_data_set(cx) + +* set up a dependent variable. + + IF (dsg_as_time) THEN + mv = mv_list( 1 ) + cx = cx_list(1) + cx_lo_ss(cx, z_dim) = cx_lo_ss(cx_list(ipl+1), t_dim) + cx_hi_ss(cx, z_dim) = cx_hi_ss(cx_list(ipl+1), t_dim) + IF (cx_hi_ss(cx, z_dim) .EQ. unspecified_int4) EXIT ! loop 200 + ENDIF + CALL EXTRACT_LINE ( cx, + . memry(mv)%ptr, + . mv, + . dep_dat, + . dep_dim, + . ndim, + . npts2, + . status ) + ENDIF + iaxis = CGRID_AXIS(dep_dim,cx) + IF (dsg_as_time) THEN + dep_dim = t_dim + cx = cx_list( ipl ) + 1 ! for indep-data, each line + mv = mv_list( 1 ) + ENDIF + + IF ( status .NE. ferr_ok ) GOTO 5000 + + IF (cxdsg_empty_set) THEN + dep_dat(1) = -1. + dep_dat(2) = 1. + ENDIF + + +* ... If this is a profile or timeseries collection maybe we are just making a map with +* the stations so dep_dat and indep dat are the longitudes and latitudes. The flag +* expand_dsg_ok false, says just plot the instance variables + +* ... Or, is this a dsg Profile dataset? Are we doing a vs plot? If so we may need to expand the +* first or second variable to the OBS axis. + + IF (its_dsg .AND. + . orientation.EQ.pfeatureType_Profile .AND. + . npts2.EQ.nfeatures .AND. expand_dsg_ok) THEN + IF (npts2 .EQ. nfeatures) CALL DSG_OBS_BY_FEATURE_VAR + . (dset, cx, nfeatures, obsdimlen, dep_dat) + npts2 = obsdimlen + dsg_fmt_grid = dsg_xlate_grid(dset) + ENDIF + +* ... For a dsg dataset, can ask to color by a feature-level variable +* such as an ID or a string-valued feature variable. + + doing_ribbon_color = ribbon_plot .AND. ipl.EQ.ribbon_var + + IF (its_dsg .AND. doing_ribbon_color) THEN + + color_dsg_id = (cx_variable(cx).EQ.dsg_feature_var(dset) ) .OR. + . (npts2.LT.obsdimlen .AND. cx_type(cx).EQ.ptype_string) + IF (color_dsg_id ) CALL SET_DSG_ID_LEV + . (dset, cx, nfeatures, dep_dat, mr_type(mv), changed_key) + + IF (npts2 .EQ. nfeatures) THEN + IF (color_dsg_id) THEN + DO i = 1, npts2 + dep_dat(i) = i + ENDDO + ENDIF + IF (expand_dsg_ok) THEN + CALL DSG_OBS_BY_FEATURE_VAR (dset, cx, nfeatures, obsdimlen, dep_dat) + npts2 = obsdimlen + ENDIF + ENDIF + + ENDIF + +* ... Is this a dsg dataset? if so apply an obs mask. Resets npts2 + + IF (grid_is_dsg .AND. .NOT.its_traj .AND. + . expand_dsg_ok .AND. .NOT.cxdsg_empty_set) + . CALL PLOT_DSG_APPLY_OBS_MASK (dset, cx, nfeatures, + . obsdimlen, dep_dat, npts2) + +* ... replicate points to create a "step" plot if requested (1/01) + nload = npts2 * step_inc + IF (step_inc .EQ. 2) THEN + DO i = npts2, 1, -1 + dep_dat(2*i) = dep_dat(i) + dep_dat(2*i-1) = dep_dat(i) + END DO + ELSEIF (step_inc .EQ. 3) THEN + DO i = npts2, 1, -1 + dep_dat(3*i) = mr_bad_data(mv) + dep_dat(3*i-1) = dep_dat(i) + dep_dat(3*i-2) = dep_dat(i) + END DO + ENDIF + +* set up corresponding independent axis +* ( note - each variable may have different points on the indep. axis ) + +* ... number of values must be equal in independ. and depend. for ordered pairs + IF ( npts2 .NE. npts ) GOTO 5130 + +* Set up the dependent axis. + + set_axis = .TRUE. + IF (ribbon_plot .OR. addgaps) set_axis = (ipl .EQ. var1) + IF (set_axis) THEN + IF (.NOT.overlay) THEN + CALL MINMAX( dep_dat, npts, mr_bad_data(mv), lo, hi, nok ) + IF (nok .EQ. 0) THEN ! all missing + dep_min = -1 + dep_max = 1 + ELSE + dep_min = MIN(dep_min, lo) + dep_max = MAX(dep_max, hi) + ENDIF + ENDIF + + all_1_dep = NO_LINE_RANGE( dep_dat, nload, mr_bad_data(mv), val1 ) + IF (.NOT.all_1_dep .AND. .NOT.dep_is_log) THEN + pad = (dep_max - dep_min)* 0.01 + dep_min = dep_min - pad + dep_max = dep_max + pad + ENDIF + + cx_y = cx_list(2) + buff1 = VAR_UNITS(cx_y) + IF (doing_ribbon_color) buff1 = ' ' + + iunits = TM_UNIT_ID(buff1) + + IF (.NOT.overlay) THEN + +* Does it have time units? If so draw a formatted time axis + y_is_time = .FALSE. + since_T0 = MAX( INDEX(buff1,'since'), INDEX(buff1,'SINCE') ) + slen = TM_LENSTR1(buff1) + IF (.NOT.x_is_time .AND. since_T0.GT.0 .AND. slen.GT.since_T0+5) THEN + cal_id_1 = 1 + CALL TM_ALLO_DYN_LINE( itmp, status ) + CALL CD_GET_TIME_UNITS ( buff1, cal_id_1, line_units(itmp), + . line_t0(itmp), dum, status ) + IF (status .EQ. ferr_ok) THEN + y_is_time = .TRUE. + time_axis = .TRUE. + the_taxis = dep_dim + has_time_axis = .TRUE. + iunits = TM_UNIT_ID( line_units(itmp) ) + line_unit_code(itmp) = TM_UNIT_ID( line_units(itmp) ) + line_tunit(itmp) = un_convert(line_unit_code(itmp)) + line_cal_name(itmp) = 'GREGORIAN' + line_direction( itmp ) = 'TI' + + grid_line(the_taxis, mgrid_buff) = itmp + + igrd_save = cx_grid(cx_list(1)) + cx_grid(cx_list(1)) = mgrid_buff + ww_save = cx_lo_ww(the_taxis,cx_list(1) ) + cx_lo_ww(the_taxis,cx_list(1)) = lo + + dep_min = TSTEP_TO_SECS (mgrid_buff, the_taxis, lo) + dep_max = TSTEP_TO_SECS (mgrid_buff, the_taxis, hi) + + IF (dsg_fmt_grid .GT. 0) THEN + dep_min = TSTEP_TO_SECS (dsg_fmt_grid, t_dim, lo) + dep_max = TSTEP_TO_SECS (dsg_fmt_grid, t_dim, hi) + ENDIF + + CALL TAXIS_STYLE( dep_ax, dep_min, dep_max, tstyle, xtra_lab ) + +* 11/16 Get the min and max as used by PPLUS. That function uses WHOI formatted time + CALL TPLOT_AXIS_ENDS (dep_min, dep_max, cal_id_1, tstyle) + + CALL AXIS_END_SYMS( dep_ax, + . SECS_TO_TSTEP( mgrid_buff, the_taxis, dep_min ), + . SECS_TO_TSTEP( mgrid_buff, the_taxis, dep_max ) ) + + ENDIF + + ENDIF + +* This routine checks the units and the setting for formatted lon/lat (time?) axes +* If the formatting has been turned off, resets iunits and flag mod_vs_x. + + CALL GEOG_LABEL_VS (buff1, iunits, y_dim, axis_y_dir) + + IF (iunits .EQ. 4) THEN + delta = unspecified_val8 + + CALL AXIS_ENDS(dep_ax,dep_dim,grid1,dep_min,dep_max, + . delta, dep_is_log, dep_axtyp, versus, status) + + CALL GET_AXIS_FORMAT( dep_min, dep_max, delta, + . fmt, use_nice ) + + IF (use_nice) THEN + +* Make sure if we are using lon/lat labeling that we have just one of each direction. +* e.g. if they said just degrees_n for the x dir, but just "degrees" for the y var, +* use east for the vert. If they said degrees_e for the y var but just degrees for the x, +* use lon for the y axis. +* If the axes are labeled with longitude, latitude, or time, then set ul_dolab to not +* put upper-left labels on the plot with the ranges. + + IF (axis_x_dir.EQ.y_dim .AND. axis_y_dir.EQ.no_dim) THEN + IF (iunits .EQ. 4) THEN + ppl_buff = 'YFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LONE'''')' + ul_dolab(x_dim) = .FALSE. + plot_axis(1) = 1 + ENDIF + + ELSEIF (axis_y_dir .EQ. x_dim) THEN + ppl_buff = 'YFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LONE'''')' + mod_vs_y = .TRUE. + ul_dolab(x_dim) = .FALSE. + plot_axis(1) = 1 + + IF (axis_x_dir .EQ. no_dim) THEN + buff1 = VAR_UNITS(cx_x) + iiunits = TM_UNIT_ID(buff1) + CALL GEOG_LABEL_VS (buff1, iiunits, x_dim, axis_x_dir) + IF (iiunits .EQ. 4) THEN + ppl_buff = 'XFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LAT'''')' + ul_dolab(y_dim) = .FALSE. + plot_axis(1) = 2 + ENDIF + ENDIF + ELSEIF (axis_y_dir .EQ. y_dim) THEN + ppl_buff = 'YFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LAT'''')' + ul_dolab(y_dim) = .FALSE. + plot_axis(2) = 2 + ELSE + ppl_buff = 'YFOR,('//fmt(:TM_LENSTR1(fmt))// + . ',''''LAT'''')' + ul_dolab(y_dim) = .FALSE. + plot_axis(2) = 2 + ENDIF + + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + ENDIF + + ELSEIF (iunits .LT. 0) THEN ! time axis formatting + CALL PPLCMD ( from, line, 0, 'TIME', 1, 1 ) + CALL PPLCMD ( from, line, 0, 'TAXIS/YAXIS 60,ON', 1, 1 ) + ppl_buff = ' ' + WRITE ( ppl_buff, 3004 ) 'GREGORIAN' + CALL PPLCMD (from, line, 0, ppl_buff, 1, 1) + ul_dolab(t_dim) = .FALSE. + plot_axis(2) = 4 + + ELSE + CALL PPLCMD ( from, line, 0, 'YFOR 0', 1, 1 ) + ENDIF + ENDIF ! not overlay + ENDIF ! set_axis + + + +* Compute the mean and standard dev for ribbon-color variable, +* needed for computing color levels. Results stored in PPLUS common. +* If histogram-based levels are requested, compute the +* histogram bins. +* send REAL*4 rbad to compare with lev_max, lev_min inside compute_mnstd +* If the user set some /LEVELS then compute_mnstd just returns + + IF (doing_ribbon_color) THEN + + IF (need_histo) THEN +* create temporary buffer to contain workspace + CALL CREATE_TEMP_MEM_VAR( cx, mvh_temp, status ) + IF ( status .NE. ferr_ok ) RETURN + plot_mem_used = plot_mem_used + npts + CALL COMPUTE_HISTO_BINS (dep_dat, memry(mvh_temp)%ptr, + . mr_bad_data(mv), npts, status) + +* ... clean up temporary variable + CALL DELETE_VARIABLE( mvh_temp ) + ELSE + rbad = bad_val4 + CALL COMPUTE_MNSTD (dep_dat, mr_bad_data(mv), need_std, npts, rbad, status) + IF (need_std .AND. status.NE.ferr_ok) THEN +* ... set up for automatic levels + CALL PPLCMD ( from, line, 0, 'LEV,()', 1, 1 ) + CALL USE_LINEAR_LEVELS + status = ferr_ok + ENDIF + + ENDIF + + ENDIF + +* ... check that there is a range of dependent data for PLOT+ auto-scaling +* ... when all the variables are considered together (but not the ribbon var) + + IF (addgaps .AND. ipl.EQ.gap_var) GOTO 199 + + IF (.NOT.ribbon_plot .OR. (ribbon_plot .AND. ipl.LT.ribbon_var)) THEN + this_no_range(ipl) = .FALSE. + IF (set_axis) this_no_range(ipl) = NO_LINE_RANGE( + . dep_dat, nload, mr_bad_data(mv), val1 ) + all_1_dep = all_1_dep .AND. this_no_range(ipl) + IF (this_no_range(ipl) .AND. val1.NE.mr_bad_data(mv) ) + . this_no_range(ipl) = .FALSE. ! Keep the value for putting NO VALID on labels + + IF ( all_1_dep ) THEN + IF ( only_val .EQ. bad_val4 ) THEN + IF ( val1 .NE. mr_bad_data(mv) ) only_val = val1 + ELSE + all_1_dep = all_1_dep .AND. only_val .EQ. val1 + ENDIF + + ENDIF + ENDIF ! checking range when not ribbon variable + + 199 CONTINUE ! checking range when not gap variable + +* pass the data to PLOT+ + + WRITE ( ppl_buff, 3005 ) mr_bad_data( mv ), dep_ax + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + IF ( time_axis ) THEN + IF ( tax .EQ. mnormal .OR. tax .EQ. munknown ) THEN ! 8/91 bug + dt_min = 1.0 + ELSE + dt_min = line_tunit( grid_line(the_taxis,grid) ) / 60.! sec-->min + IF ( ITSA_TRUEMONTH_AXIS(grid_line(the_taxis,grid)) ) + . dt_min = un_convert(pun_day)/60. + +* Check for calendar mismatch for time- overlay plots + IF (overlay) THEN + IF (cal_id_1 .NE. saved_calendar_id) THEN + cal_id_old = saved_calendar_id + cal_id_new = cal_id_1 + GO TO 5160 + ENDIF + ENDIF + t1_date = INTERNAL_WHOI_DATE( grid, the_taxis, 1.0D0 ) ! moved up one line... + ENDIF + ELSE + IF (overlay .AND. has_time_axis) THEN + dt_min = saved_dt_min + t1_date = saved_t1_date + ELSEIF (.NOT.ribbon_plot) THEN + dt_min = 1.0 + t1_date = ' ' + ENDIF + ENDIF + + IF (overlay) THEN + SOVER = .true. + ELSE + has_time_axis = time_axis + IF (has_time_axis) THEN + saved_dt_min = dt_min ! to re-use on PLOT/VS or POLYGON + saved_t1_date = t1_date + saved_calendar_id = cal_id_1 + ENDIF + ENDIF + +* put colorvar in the right place + IF (doing_ribbon_color) THEN !!! .AND. .NOT.add_ribbon + flip = .FALSE. + IF (indep_dim .EQ. z_dim) THEN + WRITE ( ppl_buff, 3005 ) mr_bad_data( mv ), indep_ax + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + ENDIF + ENDIF + +! put gap-var in the right place + IF (ipl .EQ. gap_var) THEN + flip = .FALSE. + IF (indep_dim .EQ. z_dim) THEN + WRITE ( ppl_buff, 3005 ) mr_bad_data( mv ), indep_ax + CALL PPLCMD ( from, line, 0, ppl_buff, 1, 1 ) + ENDIF + ENDIF + +* If a /vs plot has a variable with units of time, put the itmp line into +* the scratch grid mgrid_buff, so the WHOI_DATE functions can get at the +* line_* information stored above after calls to CD_GET_TIME_UNITS. + + IF (x_is_time) THEN + dt_min = line_tunit(itmp)/ 60. + + adjust_time = .FALSE. + CALL AXIS_ENDS(indep_ax,the_taxis,mgrid_buff,ind_min,ind_max, + . delta, indep_is_log, indep_axtyp, versus, status) + + tref = MERGED_WHOI_DATE( cx_list(1), the_taxis, 1, .FALSE. ) + t1_date = INTERNAL_WHOI_DATE( mgrid_buff, the_taxis, 1.0D0 ) + + cx_grid(cx_list(1)) = igrd_save + cx_lo_ww(the_taxis,cx_list(1) ) = ww_save + IF (ipl .EQ. numv) THEN + IF (itmp .NE. mnormal) CALL TM_DEALLO_DYN_LINE( itmp ) + itmp = mnormal + ENDIF + + ENDIF + + IF ( y_is_time .AND. .NOT.doing_ribbon_color) THEN + + dt_min = line_tunit(itmp)/ 60. + + adjust_time = .FALSE. + CALL AXIS_ENDS(dep_ax,the_taxis,mgrid_buff,dep_min,dep_max, + . delta, indep_is_log, indep_axtyp, versus, status) + + tref = MERGED_WHOI_DATE( cx_list(1), the_taxis, 1, .FALSE. ) + t1_date = INTERNAL_WHOI_DATE( mgrid_buff, the_taxis, 1.0D0 ) + + cx_grid(cx_list(1)) = igrd_save + cx_lo_ww(the_taxis,cx_list(1) ) = ww_save + IF (ipl .EQ. numv) THEN + IF (itmp .NE. mnormal) CALL TM_DEALLO_DYN_LINE( itmp ) + itmp = mnormal + ENDIF + + + IF (.NOT.overlay) THEN + saved_dt_min = dt_min ! to re-use on PLOT/VS or POLYGON + saved_t1_date = t1_date + saved_calendar_id = cal_id_1 + ENDIF + + IF (iunits .LT. 0) THEN + dt_min = saved_dt_min + t1_date = saved_t1_date + cal_id_1 = saved_calendar_id + ENDIF + + ENDIF + +* See pplldx comments. When loading the color-by variable lon vertical time axis +* this forces the right combination of settings. + + IF (doing_ribbon_color) icode = 2 + IF (y_is_time .AND. doing_ribbon_color) icode = 7 + + IF ( flip ) THEN + CALL PPLLDX_envelope(icode,dep_dat,indep_dat,nload, + . t1_date, tref, dt_min, plot_mem_used) + ELSE + CALL PPLLDX_envelope(icode,indep_dat,dep_dat,nload, + . t1_date, tref, dt_min, plot_mem_used) + ENDIF +* ... increment number of lines on plot + nline_on = nline_on + 1 + nline_in_mem = nline_in_mem + 1 +* ... assign line style for the data plotting: + CALL LINE_STYLE(line_symbol, sym_size, skipsym, color, color1, use_line, + . do_dash, dashstyle, nline_in_mem, nline_on) + + nload_all = MAX(nload, nload_all) + 200 CONTINUE + + nline_on = nline_on + color1 - 1 + +* axis scaling and formatting +* ... independent axis + + IF ( .NOT.overlay ) THEN + +* ... dependent axis scaling +* force axis scaling if the data has no range + + IF (the_taxis .EQ. 0) the_taxis = dep_dim + + IF ( pdeplim .GT. 0 ) THEN + IF (pylim.GT.0 .AND. .NOT.denig_xylim_msg_done) THEN + CALL WARN( '/XLIMITS and /YLIMITS are deprecated.') + CALL WARN( 'Use /HLIMITS and /VLIMITS instead.') + denig_xylim_msg_done = .TRUE. + ENDIF +* Check for valid log axis before proceeding +* First make sure the min and max are set. + IF (dep_is_log) THEN + + CALL MINMAX( dep_dat, npts, mr_bad_data(mv), lo, hi, nok ) + IF (nok .EQ. 0) GOTO 5171 + dep_min = lo + dep_max = hi + + CALL AXIS_ENDS( dep_ax, the_taxis, grid, dep_min, dep_max, + . delta, dep_is_log, dep_axtyp, versus, status ) + + ENDIF + +* Do the dependent variable(s) have units of time? If not then don't send +* dep_dim = t_dim to EQUAL_RANGE. When the dependent data is on a time +* axis, then the_taxis is set, but we don't want to be translating limits +* for the dependent axis as time units. + + dep_dir = dep_dim + IF (dep_dim .EQ. t_dim) THEN + + DO ipl = 1, nmv + buff1 = VAR_UNITS( cx_list(ipl) ) + iunits = TM_UNIT_ID(buff1) + IF ( iunits .GT. pun_last_time ) dep_dir = 0 + ENDDO + + ENDIF + + CALL EQUAL_RANGE( + . cmnd_buff(qual_start(pdeplim):qual_end(pdeplim)), + . dep_dir, dep_min, dep_max, delta, + . formatted, range_rqd, cal_id_1, status ) + IF ( status .NE. ferr_OK ) GOTO 5000 + + CALL AXIS_ENDS( dep_ax, the_taxis, grid, dep_min, dep_max, + . delta, dep_is_log, dep_axtyp, versus, status ) + + ELSEIF ( all_1_dep ) THEN + IF (only_val .EQ. bad_val4) val1 = 0.0 ! 10/99 + delta = 1. + IF (val1 .NE. 0) delta = 0.1* ABS(val1) + + CALL AXIS_ENDS( dep_ax, the_taxis, grid,val1-delta,val1+delta,delta, + . dep_is_log, dep_axtyp, versus, status ) + + ELSEIF (dep_is_log) THEN ! get dependent var scaling and set up log axis + CALL MINMAX( dep_dat, npts, mr_bad_data(mv), lo, hi, nok ) + IF (nok .EQ. 0) GOTO 5171 ! all missing + dep_min = lo + dep_max = hi + + CALL AXIS_ENDS( dep_ax, the_taxis, grid, dep_min, dep_max, + . delta, dep_is_log, dep_axtyp, versus, status ) + + ELSE + + IF (its_dsg .AND. versus) THEN ! /vs plot with depth as an axis? + dsg_fmt_grid = dsg_xlate_grid(dset) + cx = cx_list( 2 ) + iaxis = grid_line(z_dim,dsg_fmt_grid) + IF (cx_variable(cx) .EQ. dsg_coord_var(z_dim,dset) .AND. + . line_direction(iaxis) .EQ. "UD") THEN + tmp = dep_min + dep_min = dep_max + dep_max = tmp + ENDIF + + delta = unspecified_val8 + CALL AXIS_ENDS( dep_ax, the_taxis, grid, dep_min, dep_max, + . delta, dep_is_log, dep_axtyp, versus, status ) + ENDIF + + + ENDIF + IF ( status .NE. ferr_OK ) THEN + first = dep_min + last = dep_max + GOTO 5170 + ENDIF + +* When not an overlay, set the axis type to log or reverse log + + IF (dep_is_log .OR. indep_is_log) THEN + IF (flip) THEN + WRITE (val_buff, 3006) dep_axtyp, indep_axtyp + ELSE + WRITE (val_buff, 3006) indep_axtyp, dep_axtyp + ENDIF + + 3006 FORMAT ('axtype,', I2, ',', I2) + CALL PPLCMD ( from, line, 0, val_buff, 1, 1) + ENDIF + + IF (cxdsg_empty_set) THEN + WRITE ( val_buff, '(3(E14.7,1X))' ) dep_dat(1), dep_dat(2), dep_dat(2) + CALL PPLCMD ( from, line, 0, + . dep_ax//'AXIS '//val_buff , 1, 1 ) + WRITE (cxdsg_axlabp_save, "(I2, ',', I2)") LABELX, LABELY + IF (flip) THEN + CALL PPLCMD (from, line, 0,'AXLABP,0,-1',1,1) + ELSE + CALL PPLCMD (from, line, 0,'AXLABP,-1,0',1,1) + ENDIF + ENDIF + + ax1 = axdir(indep_dim) + IF (flip) THEN + CALL PPLCMD ( from, line, 0, 'SET AX_VERT '//ax1, 1, 1 ) + ELSE + CALL PPLCMD ( from, line, 0, 'SET AX_HORIZ '//ax1, 1, 1 ) + ENDIf + ENDIF + +* TITLES +* Main plot title. Add labels for multi-lines as a key, +* or a "Colored by var2" label for the ribbon color variable +* + hlen = xlen ! single-precision xlen from PPLUS common -> double prec. var + + IF (x_is_time .OR. y_is_time) THEN + ul_dolab(t_dim) = .FALSE. + READ (t1_date(13:14), *) cen + READ (t1_date(1:2), *) yr + ind_min = 100*cen + yr + ENDIF + +* Line-key labels + + IF ( .NOT.no_labels ) THEN + + CALL LINE_PLOT_LABELS (var1, nkey_entries, ndv, cx_list, + . this_no_range, overlay, versus, nokey, time_axis, + . tstyle, cal_id_1, ribbon_var, indep_lab, dep_lab, + . ind_min, dep_len, hlen, 0, dset, 0) + ENDIF + +* successful completion + 1000 CONTINUE + status = ferr_ok + IF (overlay) iautot = iautot_save + + IF (dsg_as_time) nmv = nmv + 1 ! for un-doing all the extra contexts + + RETURN + +* error exit + 5000 CALL PPLCMD ( from, line, 0, 'NLINES', 1, 1 ) ! wipe buffers clean + IF (itmp .NE. mnormal) CALL TM_DEALLO_DYN_LINE( itmp ) + IF (dsg_as_time) nmv = nmv + 1 ! for un-doing all the extra contexts + RETURN + 5050 buff1 = LEFINT( max_line_on_plot, slen ) + CALL ERRMSG( ferr_invalid_command, status, 'cannot plot more than ' + . //buff1(:slen)//' lines with a single PLOT', *5000) + 5100 CALL ERRMSG( ferr_invalid_command, status, + . cmnd_buff(:len_cmnd)//' : vs what ?', *5000 ) + 5120 CALL ERRMSG( ferr_dim_underspec, status, + . 'overlay is on a different axis'//pCR// + . '"'//cmnd_buff(:len_cmnd)//'"', *5000 ) + 5130 buff1 = LEFINT( npts, slen ) + buff2 = LEFINT( npts2, slen2 ) + buff3 = LEFINT( ipl, slen3 ) + CALL ERRMSG( ferr_dim_underspec, status, + . 'unequal line lengths: '//pCR// + . 'First expression has '//buff1(:slen)//' points.'//pCR// + . 'Expression '//buff3(:slen3)//' has '//buff2(:slen2)// + . ' points:'// + . pCR//'"'//cmnd_buff(:len_cmnd)//'"', *5000 ) + + 5140 buff3 = LEFINT( ipl, slen3 ) + CALL ERRMSG( ferr_dim_underspec, status, + . 'differing axes: '//pCR// + . 'first line is on '//ww_dim_name(indep_dim)//' axis'//pCR// + . 'line '//buff3(:slen3)//' is on '//ww_dim_name(dep_dim)// + . ' axis', *5000 ) + +! 5150 buff3 = LEFINT( INT(0.999*pplmem_nsize), slen3 ) +! buff2 = LEFINT( plot_mem_used, slen2 ) +! CALL ERRMSG( ferr_prog_limit, status, +! . 'Requested '//buff2(:slen2)//' words to plot'//pCR// +! . 'Plot buffer size is: '//buff3(:slen3), *5000 ) + + 5160 cal_name = TM_GET_CALENDAR_NAME(cal_id_old) + cal_name_new = TM_GET_CALENDAR_NAME(cal_id_new) + slen = TM_LENSTR1 (cal_name) + slen2 = TM_LENSTR1(cal_name_new) + CALL ERRMSG( ferr_inconsist_grid, status, + . 'Differing calendar axes: '//pCR// + . 'first variable is on '//cal_name(:slen)// + . ' axis'//pCR// + . 'subsequent variable is on '//cal_name_new(:slen2)// + . ' axis', *5000 ) + + 5170 buff1 = LEFT_REAL (first, '(G15.3)', slen) + buff2 = LEFT_REAL (last, '(G15.3)', slen2) + CALL ERRMSG( ferr_out_of_range, status, + . 'Limits for log axis negative or too small: '// + . buff1(:slen)// ' : '// buff2(:slen2), *5000 ) + + 5171 CALL ERRMSG( ferr_out_of_range, status, + . 'Data all-missing. Cannot define log axis for plotting', + . *5000 ) + + + 5180 CALL ERRMSG( ferr_limits, status, + . 'One-point independent axis: Requires a '// + . '/HLIMIT or /VLIMIT specification', *5000 ) + + 5190 CALL ERRMSG( ferr_syntax, status, + . 'PLOT/VS/GAPLOC only for a single dependent variable', *5000 ) + + 5200 CALL ERRMSG( ferr_grid_definition, status, + . 'Data grid is not a DSG grid ', + . *5000 ) + + 5210 CALL ERRMSG( ferr_grid_definition, status, + . 'Plot trajectoryProfile data as trajectory, '// + . 'only for a single variable', + . *5000 ) + END diff --git a/fer/utl/make_dsg_feature_mask_dirs.F b/fer/utl/make_dsg_feature_mask_dirs.F new file mode 100644 index 000000000..4592c4f23 --- /dev/null +++ b/fer/utl/make_dsg_feature_mask_dirs.F @@ -0,0 +1,231 @@ + SUBROUTINE MAKE_DSG_FEATURE_MASK_DIRS (dset, cx, fmask, nfeatures, dirs ) + +* +* +* This software was developed by the Thermal Modeling and Analysis +* Project(TMAP) of the National Oceanographic and Atmospheric +* Administration's (NOAA) Pacific Marine Environmental Lab(PMEL), +* hereafter referred to as NOAA/PMEL/TMAP. +* +* Access and use of this software shall impose the following +* obligations and understandings on the user. The user is granted the +* right, without any fee or cost, to use, copy, modify, alter, enhance +* and distribute this software, and any derivative works thereof, and +* its supporting documentation for any purpose whatsoever, provided +* that this entire notice appears in all copies of the software, +* derivative works and supporting documentation. Further, the user +* agrees to credit NOAA/PMEL/TMAP in any publications that result from +* the use of this software or in any product that includes this +* software. The names TMAP, NOAA and/or PMEL, however, may not be used +* in any advertising or publicity to endorse or promote any products +* or commercial entity unless specific written permission is obtained +* from NOAA/PMEL/TMAP. The user also understands that NOAA/PMEL/TMAP +* is not obligated to provide the user with any support, consulting, +* training or assistance of any kind with regard to the use, operation +* and performance of this software nor to provide the user with any +* updates, revisions, new versions or "bug fixes". +* +* THIS SOFTWARE IS PROVIDED BY NOAA/PMEL/TMAP "AS IS" AND ANY EXPRESS +* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +* ARE DISCLAIMED. IN NO EVENT SHALL NOAA/PMEL/TMAP BE LIABLE FOR ANY SPECIAL, +* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER +* RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF +* CONTRACT, NEGLIGENCE OR OTHER TORTUOUS ACTION, ARISING OUT OF OR IN +* CONNECTION WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE. +* +* +* Create a feature mask for a DSG based upon world coordinate limits +* in the given context and anhy masks that are set. The argument dirs(4) +* lets set which directions to use in applying the context + +* NOAA/PMEL, Seattle, WA - Tropical Modeling and Analysis Program + +* V720 4/17 *sh* +* V73+ *acm* 12/17 INIT_CXDSG needs nfeatures +* V741 6/18 *acm* Mask check enhanced for partial overlap of feature and region +* and handle modulo longitudes +* V741 7/18 *acm* Fixes for modulo longitudes in computing whether context +* and coordinate regions overlap +* V74+ 8/18 *acm* More on modulo longitudes and the mask +* V76 1/20 *acm* working with Point-type dsg data + + include 'tmap_dims.parm' + include 'ferret.parm' + include 'errmsg.parm' + include 'xtm_grid.cmn_text' + include 'xdset_info.cmn_text' + include 'xdyn_linemem.cmn_text' + include 'xcontext.cmn' + include 'xdsg_context.cmn' + include 'xprog_state.cmn' + include 'xvariables.cmn' +# include "tmap_dset.parm" + + +* calling argument declarations: + INTEGER dset, cx, nfeatures + LOGICAL fmask(nfeatures), dirs(*) + +* internal variable declarations: + LOGICAL one_or_other, ismod, its_dsg, its_cmpnd, ignore_depth + INTEGER TM_DSG_NF2FEATURES, + . grid, i, idim, lm, orientation, obsdimlen, feature_line, + . nftrsets, n2, lm_index + REAL GET_LINE_COORD, coord, ftr_lo(4,nfeatures), ftr_hi(4,nfeatures), + . flo, fhi, axmod, bad, val + + + INTEGER, DIMENSION(:), ALLOCATABLE :: station_index + LOGICAL, DIMENSION(:), ALLOCATABLE :: f2mask + + + axmod = 360.0 + +* Will allow depth coordinates to be missing for timeseries or trajectory data + + grid = dsg_xlate_grid(dset) + CALL TM_DSG_FACTS( grid, orientation, obsdimlen, feature_line, + . its_dsg, its_cmpnd ) + ignore_depth = orientation.EQ.pfeatureType_Trajectory .OR. + . orientation.EQ.pfeatureType_Point .OR. + . orientation.EQ.pfeatureType_TimeSeries + +! DSGTBD: presently the same constraints mask gets re-computed many times: +! ... during each IS_DO_OP, during regridding, and for final LIST, PLOT, etc. +! It would be desirable to avoid this. A way to do this would be +! 1. compute a unique hash based upon the XYZTE context constraints +! 2. save the masks that computed here in lm line memory +! 3. use the hash value and dataset number to recover the saved masks +! 4. clean out the stored masks (say?) when the dataset gets closed + +* initialize + cxdsg_constrain_e = .FALSE. + IF (cx .GT. cx_none) CALL INIT_CXDSG(dset, cx, nfeatures) ! summarize context ww limits + +* preset the features and obs masks based upon E constraints + IF (cxdsg_constrain_e) THEN + DO i = 1, nfeatures + fmask(i) = i .GE. cxdsg_lo_e .AND. i .LE. cxdsg_hi_e + ENDDO + ELSE + DO i = 1, nfeatures + fmask(i) = .TRUE. + ENDDO + ENDIF + +* Did the user send in a station or trajectory mask with SET DATA/TRAJMASK= or SET DATA/STNMASK= + + IF (dsg_ftrsetmsk_lm(dset).NE.unspecified_int4 .AND. + . dsg_ftrsetmsk_lm(dset).NE.int4_init ) THEN + + nftrsets = TM_DSG_NF2FEATURES (dset) + ALLOCATE (f2mask(nftrsets)) + + ALLOCATE (station_index(nfeatures)) + + DO i = 1, nftrsets + val = GET_LINE_COORD (linemem(dsg_ftrsetmsk_lm(dset))%ptr, i) + f2mask(i) = (val.NE.0. .AND. val.NE.bad_val4) + ENDDO + + lm_index = dsg_loaded_lm(dsg_index_var(dset)) + DO i = 1, nfeatures + station_index(i) = dsg_linemem(lm_index)%ptr(i) + 1 + ENDDO + +* which feature is in the station or traj + + DO i = 1, nfeatures + n2 = station_index(i) + fmask(i) = f2mask(n2) + ENDDO + + DEALLOCATE(f2mask) + DEALLOCATE(station_index) + +* Did the user send in a mask with e.g. "PLOT/FMASK var, maskvals" + + ELSE IF (dsg_msk_lm(dset) .NE. unspecified_int4) THEN + DO i = 1, nfeatures + val = GET_LINE_COORD (linemem(dsg_msk_lm(dset))%ptr, i) + fmask(i) = (val.NE.0. .AND. val.NE.bad_val4) + ENDDO + ENDIF + +* quick exit? + IF (cxdsg_no_coord_constraints .OR. cx .EQ. cx_none) RETURN + +* get the coordinate extremes feature by feature + CALL DSG_FEATURE_LIMS(dset, nfeatures, ftr_lo, ftr_hi) + +* first mask at the feature level based on instance coordinates + DO idim = 1, t_dim + IF (.NOT.cxdsg_constrain(idim)) CYCLE + IF (.NOT.cxdsg_has_coord(idim)) CYCLE + IF (.NOT.dirs(idim)) CYCLE + + bad = cxdsg_bad_val(idim) + IF (cxdsg_is_obs_coord(idim)) THEN + +* ... use the coordinate extremes for this feature + + DO i = 1, nfeatures + + IF (.NOT.fmask(i)) CYCLE + + flo = ftr_lo(idim,i) + fhi = ftr_hi(idim,i) + + IF ((flo.EQ.unspecified_val8 .OR. fhi.EQ.unspecified_val8) .AND. + . (idim.EQ.z_dim.OR.idim.EQ.t_dim) .AND. + . ignore_depth ) CYCLE ! Depths might be all bad, feature still ok + + IF (idim .EQ. x_dim) THEN + CALL MODSCAT ( cxdsg_constrain_lo(idim), + . cxdsg_constrain_hi(idim), axmod, 1, flo) + CALL MODSCAT ( cxdsg_constrain_lo(idim), + . cxdsg_constrain_hi(idim), axmod, 1, fhi) + ismod = flo.NE.ftr_lo(idim,i) .OR. fhi.NE.ftr_hi(idim,i) + ENDIF + +* Does the requested region lie partially or entirely within the feature? + one_or_other = + . ( (cxdsg_constrain_lo(idim) .GE. flo .AND. + . cxdsg_constrain_lo(idim) .LE. fhi ) .OR. + . (cxdsg_constrain_hi(idim) .GE. flo .AND. + . cxdsg_constrain_hi(idim) .LE. fhi) ) + +* or if the feature lies within the requested region that's ok too + IF (.NOT.one_or_other) THEN + one_or_other = + . ( (flo .GE. cxdsg_constrain_lo(idim) .AND. + . flo .LE. cxdsg_constrain_hi(idim) ) .OR. + . (fhi .GE. cxdsg_constrain_lo(idim) .AND. + . fhi .LE. cxdsg_constrain_hi(idim)) ) + ENDIF + + fmask(i) = fmask(i).AND. one_or_other + ENDDO + + ELSE + +* ... use the single point coordinate for this feature + + lm = cxdsg_coord_lm(idim) + DO i = 1, nfeatures + coord = dsg_linemem(lm)%ptr(i) + IF (coord .EQ. bad) CYCLE + IF (idim .EQ. x_dim) CALL MODSCAT ( + . cxdsg_constrain_lo(idim), + . cxdsg_constrain_hi(idim), axmod, 1, coord) + + fmask(i) = fmask(i) + . .AND. coord .LE. cxdsg_constrain_hi(idim) + . .AND. coord .GE. cxdsg_constrain_lo(idim) + ENDDO + ENDIF + ENDDO + + RETURN + END diff --git a/fer/utl/make_dsg_ftrset_mask.F b/fer/utl/make_dsg_ftrset_mask.F new file mode 100644 index 000000000..2fdb6199d --- /dev/null +++ b/fer/utl/make_dsg_ftrset_mask.F @@ -0,0 +1,98 @@ + SUBROUTINE MAKE_DSG_FTRSET_MASK(dset, cx, nfeatures, fmask, + . nftrsets, smask ) +* +* +* This software was developed by the Thermal Modeling and Analysis +* Project(TMAP) of the National Oceanographic and Atmospheric +* Administration's (NOAA) Pacific Marine Environmental Lab(PMEL), +* hereafter referred to as NOAA/PMEL/TMAP. +* +* Access and use of this software shall impose the following +* obligations and understandings on the user. The user is granted the +* right, without any fee or cost, to use, copy, modify, alter, enhance +* and distribute this software, and any derivative works thereof, and +* its supporting documentation for any purpose whatsoever, provided +* that this entire notice appears in all copies of the software, +* derivative works and supporting documentation. Further, the user +* agrees to credit NOAA/PMEL/TMAP in any publications that result from +* the use of this software or in any product that includes this +* software. The names TMAP, NOAA and/or PMEL, however, may not be used +* in any advertising or publicity to endorse or promote any products +* or commercial entity unless specific written permission is obtained +* from NOAA/PMEL/TMAP. The user also understands that NOAA/PMEL/TMAP +* is not obligated to provide the user with any support, consulting, +* training or assistance of any kind with regard to the use, operation +* and performance of this software nor to provide the user with any +* updates, revisions, new versions or "bug fixes". +* +* THIS SOFTWARE IS PROVIDED BY NOAA/PMEL/TMAP "AS IS" AND ANY EXPRESS +* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +* ARE DISCLAIMED. IN NO EVENT SHALL NOAA/PMEL/TMAP BE LIABLE FOR ANY SPECIAL, +* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER +* RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF +* CONTRACT, NEGLIGENCE OR OTHER TORTUOUS ACTION, ARISING OUT OF OR IN +* CONNECTION WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE. +* +* +* Create a timeseries-station or trajectory mask for a DSG based upon +* world coordinate limits and/or a user-defined mask. Uses a feature-mask +* defined by a call to MAKE_DSG_FEATURE_MASK. + +* NOAA/PMEL, Seattle, WA - SDIG + +* V762 8/20 *ACM* +* V762 9/20 *ACM* for both TimeseriesProfile and TrajectoryProfile data types + + include 'tmap_dims.parm' + include 'ferret.parm' + include 'xtm_grid.cmn_text' + include 'xdset_info.cmn_text' + include 'xdyn_linemem.cmn_text' + +#include "tmap_dset.parm" + +* calling argument declarations: + INTEGER dset, cx, nfeatures, nftrsets + LOGICAL fmask(nfeatures), smask(*) + +* internal variable declarations: + LOGICAL its_dsg, its_cmpnd + INTEGER grid, i, lm, orientation, obsdimlen, feature_line, + . lm_index + + INTEGER, DIMENSION(:), ALLOCATABLE :: station_index + +* Initialize + DO i = 1, nftrsets + smask(i) = .FALSE. + ENDDO + +* Get dataset info + + grid = dsg_xlate_grid(dset) + CALL TM_DSG_FACTS( grid, orientation, obsdimlen, feature_line, + . its_dsg, its_cmpnd ) + + IF (.NOT. its_cmpnd) RETURN + +* Get the station-index and use it to pick out the station-id or trajectory-id mask +* from the feature mask + + ALLOCATE (station_index(nfeatures)) + lm_index = dsg_loaded_lm(dsg_index_var(dset)) + + DO i = 1, nfeatures + station_index(i) = dsg_linemem(lm_index)%ptr(i) + 1 + ENDDO + +* which station/traj are represented in the feature mask + + DO i = 1, nfeatures + IF (fmask(i)) smask(station_index(i)) = .TRUE. + ENDDO + + DEALLOCATE(station_index) + + RETURN + END