* * eof_stat.F * * 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. * * Ansley Manke * July 1999 * Aug 2000 change fcn description; argument can be fcn of Z as well as x,y,t * May 2001 change to simpler solver EOFIN (time series all need to be filled * * Jul 2001 Move to statically linked code in fer/efi. * Jul 2001 move count_eof and pack_eof to file eofsubs.f, rename solve_eof * to solve_eof_stat. * Make the function work w/ data having multiple depths: compute * EOF solution for each depth. * Aug 2001 move to fer/efi directory, to be statically linked to Ferret. * Change INCLUDE statements to remove directory spec. * Compute EOF/s and time amplitude funcions from a 2-d field. * Return statistics: number of eof's computed, percent variance explained * by each EOF, and the eigenvalue for each EOF. * x-axis is abstract length nx*ny; listing the statistics for each EOF. * y-axis is length 3: * 1st is the # of eof's computed * 2nd is the Pct variation for each eof * 3rd is the eigenvalue for each eof. * * Compute EOF/s and time amplitude funcions from a 2-d field. Based on * program (coadseof.f, etc.) by Dai McClurg and Ansley Manke and calls * Billy Kessler's method for finding EOFs * of NON-gappy time series EOFIN.for * * * In this SUBROUTINE we provide information about * the function. The user configurable information * consists of the following: * * descr Text description of the function * * num_args Required number of arguments * * axis_inheritance Type of axis for the result * ( CUSTOM, IMPLIED_BY_ARGS, NORMAL, ABSTRACT ) * CUSTOM - user defined axis * IMPLIED_BY_ARGS - same axis as the incoming argument * NORMAL - the result is normal to this axis * ABSTRACT - an axis which only has index values * * piecemeal_ok For memory optimization: * axes where calculation may be performed piecemeal * ( YES, NO ) * * * For each argument we provide the following information: * * name Text name for an argument * * unit Text units for an argument * * desc Text description of an argument * * axis_influence Are this argument's axes the same as the result grid? * ( YES, NO ) * * axis_extend How much does Ferret need to extend arg limits relative to result * SUBROUTINE eof_stat_init(id) INCLUDE 'ferret_cmn/EF_Util.cmn' INTEGER id, arg ************************************************************************ * USER CONFIGURABLE PORTION | * | * V CHARACTER*100 fcn_desc 10 FORMAT ('EOF statistics from XYT field. j=1:#EOFs, ', . 'j=2:%variation, j=3:eigenvalues' ) WRITE (fcn_desc, 10) CALL ef_set_desc(id, fcn_desc) CALL ef_set_num_args(id, 2) CALL ef_set_axis_inheritance(id, ABSTRACT, ABSTRACT, . IMPLIED_BY_ARGS, NORMAL) CALL ef_set_piecemeal_ok(id, NO, NO, NO, NO) CALL ef_set_num_work_arrays(id, 9) arg = 1 CALL ef_set_arg_name(id, arg, 'A') CALL ef_set_arg_unit(id, arg, ' ') CALL ef_set_arg_desc(id, arg,'Variable in x,y,t; may be fcn of z') CALL ef_set_axis_influence(id, arg, YES, YES, YES, YES) arg = 2 CALL ef_set_arg_name(id, arg, 'frac_timeser') CALL ef_set_arg_unit(id, arg, ' ') CALL ef_set_arg_desc(id, arg, . 'Use only those time series with this fraction valid data.') CALL ef_set_axis_influence(id, arg, NO, NO, NO, NO) * ^ * | * USER CONFIGURABLE PORTION | ************************************************************************ RETURN END * * In this subroutine we provide information about the lo and hi * limits associated with each abstract or custom axis. The user * configurable information consists of the following: * * loss lo subscript for an axis * * hiss hi subscript for an axis * SUBROUTINE eof_stat_result_limits(id) INCLUDE 'ferret_cmn/EF_Util.cmn' INCLUDE 'ferret_cmn/EF_mem_subsc.cmn' INTEGER id INTEGER arg_lo_ss(4,EF_MAX_ARGS), arg_hi_ss(4,EF_MAX_ARGS), . arg_incr(4,EF_MAX_ARGS) * ********************************************************************** * USER CONFIGURABLE PORTION | * | * V INTEGER my_lo_l, my_hi_l INTEGER nx, ny * Use utility functions to get context information about the * 1st argument, to set the abstract axis lo and hi indices. CALL ef_get_arg_subscripts(id, arg_lo_ss, arg_hi_ss, arg_incr) nx = arg_hi_ss(X_AXIS, ARG1) - arg_lo_ss(X_AXIS, ARG1) + 1 ny = arg_hi_ss(Y_AXIS, ARG1) - arg_lo_ss(Y_AXIS, ARG1) + 1 my_lo_l = 1 my_hi_l = nx * ny CALL ef_set_axis_limits(id, X_AXIS, my_lo_l, my_hi_l) CALL ef_set_axis_limits(id, Y_AXIS, 1, 3) * ^ * | * USER CONFIGURABLE PORTION | * ********************************************************************** RETURN END * * In this subroutine we request an amount of storage to be supplied * by Ferret and passed as an additional argument. * SUBROUTINE eof_stat_work_size(id) INCLUDE 'ferret_cmn/EF_Util.cmn' INCLUDE 'ferret_cmn/EF_mem_subsc.cmn' INTEGER id * ********************************************************************** * USER CONFIGURABLE PORTION | * | * V * * Set the work arrays, X/Y/Z/T dimensions * * ef_set_work_array_dims(id,array #,xlo,ylo,zlo,tlo,xhi,yhi,zhi,thi) * COMMON /eof_statSTOR/ mx INTEGER mx1, my1, mt1, mx, mxmt, mxmx INTEGER iwork INTEGER arg_lo_ss(4,1:EF_MAX_ARGS), arg_hi_ss(4,1:EF_MAX_ARGS), . arg_incr(4,1:EF_MAX_ARGS) CALL ef_get_arg_subscripts(id, arg_lo_ss, arg_hi_ss, arg_incr) mx1 = 1 + arg_hi_ss(X_AXIS,ARG1) - arg_lo_ss(X_AXIS,ARG1) my1 = 1 + arg_hi_ss(Y_AXIS,ARG1) - arg_lo_ss(Y_AXIS,ARG1) mt1 = 1 + arg_hi_ss(T_AXIS,ARG1) - arg_lo_ss(T_AXIS,ARG1) mx = mx1 * my1 + 10 mxmt = mx * mt1 + 10 mxmx = mx* mx + 10 * val iwork = 1 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mx, 1, 1, 1) * taf iwork = 2 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mxmt, 1, 1, 1) * pct iwork = 3 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mx, 1, 1, 1) * vec iwork = 4 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mxmx, 1, 1, 1) * c iwork = 5 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mxmx, 1, 1, 1) * ddat_1d iwork = 6 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mxmt, 1, 1, 1) * isave_jsave iwork = 7 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mx, 2, 1, 1) * ok iwork = 8 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mx1, my1, 1, 1) * eofwork iwork = 9 CALL ef_set_work_array_dims (id, iwork, 1, 1, 1, 1, . mx, 4, 1, 1) * ^ * | * USER CONFIGURABLE PORTION | * ********************************************************************** RETURN END * * In this SUBROUTINE we compute the result * SUBROUTINE eof_stat_compute(id, arg_1, arg_2, result, . val, taf, pct, vec, c, ddat_1d, isave_jsave, ok, eofwork) * arg_1 variable, function of (x,y,t) * result(i, 1,. ,. ) # eigenfunctions. * result(i, 2,. ,. ) pct variation for eigenfuncton i. * result(i, 3,. ,. ) eigenvalue for eigenfuncton i. * Work arrays: * val, taf, pct, vec, c, ddat_1d, isave_jsave, ok * val(neof) eigenvalues (Lambda) * taf(neof,nt) time amplitude functions (V). Dimensionless. * pct(neof) % variance represented by each EOF. * c(neof,neof) work space for cov matrix (garbage output) * isave_jsave save the locations of the data in the x-y plane * eofwork used by QRSYM and other solver routines; replaces * original ALPHA, BETA, BB, and P arrays. INCLUDE 'ferret_cmn/EF_Util.cmn' INCLUDE 'ferret_cmn/EF_mem_subsc.cmn' INTEGER id COMMON /eof_statSTOR/ mx INTEGER mx REAL bad_flag(EF_MAX_ARGS), bad_flag_result REAL arg_1(mem1lox:mem1hix, mem1loy:mem1hiy, mem1loz:mem1hiz, . mem1lot:mem1hit) REAL arg_2(mem2lox:mem2hix, mem2loy:mem2hiy, mem2loz:mem2hiz, . mem2lot:mem2hit) REAL result(memreslox:memreshix, memresloy:memreshiy, . memresloz:memreshiz, memreslot:memreshit) * Dimension the work arrays. REAL val(wrk1lox:wrk1hix, wrk1loy:wrk1hiy, . wrk1loz:wrk1hiz, wrk1lot:wrk1hit) REAL taf(wrk2lox:wrk2hix, wrk2loy:wrk2hiy, . wrk2loz:wrk2hiz, wrk2lot:wrk2hit) REAL pct(wrk3lox:wrk3hix, wrk3loy:wrk3hiy, . wrk3loz:wrk3hiz, wrk3lot:wrk3hit) REAL vec(wrk4lox:wrk4hix, wrk4loy:wrk4hiy, . wrk4loz:wrk4hiz, wrk4lot:wrk4hit) REAL c(wrk5lox:wrk5hix, wrk5loy:wrk5hiy, . wrk5loz:wrk5hiz, wrk5lot:wrk5hit) REAL ddat_1d(wrk6lox:wrk6hix, wrk6loy:wrk6hiy, . wrk6loz:wrk6hiz, wrk6lot:wrk6hit) REAL isave_jsave(wrk7lox:wrk7hix, wrk7loy:wrk7hiy, . wrk7loz:wrk7hiz, wrk7lot:wrk7hit) REAL ok(wrk8lox:wrk8hix, wrk8loy:wrk8hiy, . wrk8loz:wrk8hiz, wrk8lot:wrk8hit) REAL eofwork(wrk9lox:wrk9hix, wrk9loy:wrk9hiy, . wrk9loz:wrk9hiz, wrk9lot:wrk9hit) * After initialization, the 'res_' arrays contain indexing information * for the result axes. The 'arg_' arrays will contain the indexing * information for each variable's axes. INTEGER res_lo_ss(4), res_hi_ss(4), res_incr(4) INTEGER arg_lo_ss(4,EF_MAX_ARGS), arg_hi_ss(4,EF_MAX_ARGS), . arg_incr(4,EF_MAX_ARGS) ************************************************************************ * USER CONFIGURABLE PORTION | * | * V INTEGER neof, ier CHARACTER*255 err_msg, err_out INTEGER nx, ny, nt, i, j, k, l, k1 REAL frac_timeser * Initialize work arrays. DO 701 i = wrk1lox,wrk1hix DO 701 j = wrk1loy,wrk1hiy DO 701 k = wrk1loz,wrk1hiz DO 701 l = wrk1lot,wrk1hit val(i,j,k,l) = 0. 701 CONTINUE DO 702 i = wrk2lox,wrk2hix DO 702 j = wrk2loy,wrk2hiy DO 702 k = wrk2loz,wrk2hiz DO 702 l = wrk2lot,wrk2hit taf(i,j,k,l) = 0. 702 CONTINUE DO 703 i = wrk3lox,wrk3hix DO 703 j = wrk3loy,wrk3hiy DO 703 k = wrk3loz,wrk3hiz DO 703 l = wrk3lot,wrk3hit pct(i,j,k,l) = 0. 703 CONTINUE DO 704 i = wrk4lox,wrk4hix DO 704 j = wrk4loy,wrk4hiy DO 704 k = wrk4loz,wrk4hiz DO 704 l = wrk4lot,wrk4hit vec(i,j,k,l) = 0. 704 CONTINUE DO 705 i = wrk5lox,wrk5hix DO 705 j = wrk5loy,wrk5hiy DO 705 k = wrk5loz,wrk5hiz DO 705 l = wrk5lot,wrk5hit c(i,j,k,l) = 0. 705 CONTINUE DO 706 i = wrk6lox,wrk6hix DO 706 j = wrk6loy,wrk6hiy DO 706 k = wrk6loz,wrk6hiz DO 706 l = wrk6lot,wrk6hit ddat_1d(i,j,k,l) = 0. 706 CONTINUE DO 707 i = wrk7lox,wrk7hix DO 707 j = wrk7loy,wrk7hiy DO 707 k = wrk7loz,wrk7hiz DO 707 l = wrk7lot,wrk7hit isave_jsave(i,j,k,l) = 0. 707 CONTINUE DO 708 i = wrk8lox,wrk8hix DO 708 j = wrk8loy,wrk8hiy DO 708 k = wrk8loz,wrk8hiz DO 708 l = wrk8lot,wrk8hit ok(i,j,k,l) = 0. 708 CONTINUE DO 709 i = wrk9lox,wrk9hix DO 709 j = wrk9loy,wrk9hiy DO 709 k = wrk9loz,wrk9hiz DO 709 l = wrk9lot,wrk9hit eofwork(i,j,k,l) = 0. 709 CONTINUE CALL ef_get_res_subscripts(id, res_lo_ss, res_hi_ss, res_incr) CALL ef_get_arg_subscripts(id, arg_lo_ss, arg_hi_ss, arg_incr) CALL ef_get_bad_flags(id, bad_flag, bad_flag_result) * get array sizes nx = 1 + arg_hi_ss(X_AXIS,ARG1) - arg_lo_ss(X_AXIS,ARG1) ny = 1 + arg_hi_ss(Y_AXIS,ARG1) - arg_lo_ss(Y_AXIS,ARG1) nt = 1 + arg_hi_ss(T_AXIS,ARG1) - arg_lo_ss(T_AXIS,ARG1) * Get time percent parameter. frac_timeser = arg_2(arg_lo_ss(X_AXIS,ARG2), . arg_lo_ss(Y_AXIS,ARG2), arg_lo_ss(Z_AXIS,ARG2), . arg_lo_ss(T_AXIS,ARG2)) * Compute EOF for each depth k1 = arg_lo_ss(Z_AXIS,ARG1) DO 100 k = res_lo_ss(Z_AXIS), res_hi_ss(Z_AXIS) * Find the number of eofs to solve for. * neof = number of (x,y) points with frac_timeser good data. * Set OK to mark where they are. (note this also allows for * working around continental boundaries or other areas where * entire time series are missing) CALL count_neof (id, arg_1, neof, ok, nx, ny, k1, nt, mx, . bad_flag(ARG1), frac_timeser, err_msg, ier) IF (ier .ne. 0) then GOTO 5010 ENDIF * Put the data into the array ddat_1d(neof,nt) CALL pack_ef (id, arg_1, ddat_1d, isave_jsave, neof, . ok, frac_timeser, nx, ny, k1, nt) * Solve for the EOF's: eigenvectors, time functions, percent variance explained. CALL solve_eof_stat (ddat_1d, neof, nt, k, val, vec, . taf, pct, c, result, isave_jsave, eofwork, bad_flag(ARG1), . bad_flag_result, frac_timeser, err_msg, ier) IF (ier .ne. 0) GOTO 5020 k1 = k1 + arg_incr(Z_AXIS,ARG1) 100 CONTINUE RETURN 5010 CALL ef_bail_out (id, err_msg) RETURN 5020 WRITE (err_out, 5500) k1, err_msg 5500 FORMAT ('k=', I3, A250) CALL ef_bail_out (id, err_out) RETURN END SUBROUTINE solve_eof_stat (ddat_1d, neof, nt, k, val, vec, . taf, pct, c, result, isave_jsave, eofwork, bad_flag_dat, . bad_flag_result, frac_timeser, err_msg, ier) INCLUDE 'ferret_cmn/EF_Util.cmn' INCLUDE 'ferret_cmn/EF_mem_subsc.cmn' INTEGER id INTEGER res_lo_ss(4), res_hi_ss(4), res_incr(4) REAL result(memreslox:memreshix, memresloy:memreshiy, . memresloz:memreshiz, memreslot:memreshit) REAL bad_flag_dat, bad_flag_result, frac_timeser INTEGER neof, nt, ier INTEGER i, j, k, l REAL isave_jsave(wrk7lox:wrk7hix, wrk7loy:wrk7hiy, . wrk7loz:wrk7hiz, wrk7lot:wrk7hit) REAL eofwork(wrk9lox:wrk9hix, wrk9loy:wrk9hiy, . wrk9loz:wrk9hiz, wrk9lot:wrk9hit) REAL val(*) REAL taf(neof, nt) REAL pct(*) REAL vec(neof, neof) REAL c(neof, neof) REAL ddat_1d(neof,nt) CHARACTER*(*) err_msg CALL ef_get_res_subscripts(id, res_lo_ss, res_hi_ss, res_incr) IF (frac_timeser .LT. 1.) THEN CALL EOFIN_CHEL_GAP (ddat_1d, neof, nt, val, vec, taf, pct, c, . eofwork, bad_flag_dat, bad_flag_result, err_msg, ier) ELSE CALL EOFIN(ddat_1d,neof,NT,VAL,VEC,TAF,PCT,C, . eofwork, bad_flag_dat, bad_flag_result) ENDIF IF (ier .ne. 0) GOTO 999 * Result for j = 1 is the # eigenfunctions. j = res_lo_ss(Y_AXIS) DO 130 i = res_lo_ss(X_AXIS), res_hi_ss(X_AXIS) DO 110 l = res_lo_ss(T_AXIS), res_hi_ss(T_AXIS) result(i,j,k,l) = neof 110 CONTINUE 130 CONTINUE * Result for j = 2 is the percent variation explained by the eigenfunctions j = j + res_incr(Y_AXIS) DO 230 i = res_lo_ss(X_AXIS), res_hi_ss(X_AXIS) DO 210 l = res_lo_ss(T_AXIS), res_hi_ss(T_AXIS) result(i,j,k,l) = pct(i) 210 CONTINUE 230 CONTINUE * Result for j = 3 is the eigenvalue for each eigenfunction j = j + res_incr(Y_AXIS) DO 330 i = res_lo_ss(X_AXIS), res_hi_ss(X_AXIS) DO 310 l = res_lo_ss(T_AXIS), res_hi_ss(T_AXIS) result(i,j,k,l) = val(i) 310 CONTINUE 330 CONTINUE 999 CONTINUE RETURN END