From 5d68dcb2880a51b365d746c63462e29c803b9e1e Mon Sep 17 00:00:00 2001 From: Lawrence Takacs Date: Wed, 5 Jun 2024 10:39:28 -0400 Subject: [PATCH 1/3] Updates for automated TEM Diagnostics --- plots/grads_util/Create_TEM_ddf.csh | 95 ++ plots/grads_util/getinfo | 16 + plots/grads_util/getobs.gs | 1 + plots/grads_util/makezdif.gs | 4 +- plots/grads_util/rgbset.gs | 9 + plots/grads_util/writegrads.gs | 105 +- plots/hcmp/hcmp2 | 15 +- plots/portrait.script | 5 +- plots/quickplot | 13 +- plots/res/VERIFICATION.MERRA-2.rc | 4 +- plots/res/VERIFICATION.rc.tmpl | 4 +- plots/res/epflx.gs | 48 +- plots/res/epflx_diff.gs | 311 ++++-- plots/res/plot_season.gs | 63 +- plots/res/res | 10 + plots/res/res_2.gs | 60 +- plots/res/res_3.gs | 72 +- plots/res/setup_epflx.gs | 187 +++- plots/res/setup_wstar.gs | 25 +- plots/res/turn_around_lats.gs | 77 +- plots/res/zonal.gs | 400 ++++--- plots/zcmp/zcloseplt | 16 +- plots/zcmp/zplt | 8 +- post/CMakeLists.txt | 6 +- post/Create_TEM_Diag | 102 ++ post/Create_TEM_Diag.F | 212 ++++ post/TEM.F90 | 1575 +++++++++++++++++++++++++++ post/gcmclim.script | 4 +- post/gcmpost.script | 178 ++- post/time_ave.F | 132 ++- 30 files changed, 3345 insertions(+), 412 deletions(-) create mode 100755 plots/grads_util/Create_TEM_ddf.csh create mode 100755 post/Create_TEM_Diag create mode 100644 post/Create_TEM_Diag.F create mode 100644 post/TEM.F90 diff --git a/plots/grads_util/Create_TEM_ddf.csh b/plots/grads_util/Create_TEM_ddf.csh new file mode 100755 index 00000000..c588923a --- /dev/null +++ b/plots/grads_util/Create_TEM_ddf.csh @@ -0,0 +1,95 @@ +#!/bin/csh -fx + +setenv ARCH `uname` + +setenv SITE NCCS +setenv GEOSDIR /discover/nobackup/ltakacs/TAGS/Jason-4_0_p2/GEOSagcm +setenv GEOSBIN /discover/nobackup/ltakacs/TAGS/Jason-4_0_p2/GEOSagcm/Linux/bin +setenv GEOSUTIL /discover/nobackup/ltakacs/TAGS/Jason-4_0_p2/GEOSagcm/src/GMAO_Shared/GEOS_Util + +source $GEOSBIN/g5_modules +setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:${BASEDIR}/${ARCH}/lib + +set EXPID = $1 +set TEM_Collection = $2 + +# Count Number of "Dots" in EXPID +# ------------------------------- +set expdots = '' + @ n = 0 + @ b = 1 + set bit = `echo ${EXPID} | cut -b $b` + while( "${bit}" != '' ) + if( "${bit}" == '.' ) then + @ n = $n + 1 + endif + @ b = $b + 1 + set bit = `echo ${EXPID} | cut -b $b` + end + set expdots = `echo $expdots $n` + @ date_node = $expdots + 4 +echo Date_Node: $date_node + +/bin/rm -f grads.commands +touch grads.commands +echo \'open residual.$EXPID.ctl\' >> grads.commands +echo \'run writegrads.gs -name $EXPID.${TEM_Collection}.monthly -time one\' >> grads.commands +echo \'quit\' >> grads.commands +grads -b -l -c "run grads.commands" + +set files = `/bin/ls -1 $EXPID.${TEM_Collection}.monthly.*.data` +foreach file ($files) +set date = `echo $file | cut -d. -f${date_node}` +$GEOSBIN/flat2hdf.x -flat $file -ctl $EXPID.${TEM_Collection}.monthly.$date.ctl -nymd $date -nhms 0 -ndt 21600 +end + +stripname data.nc4 nc4 +/bin/rm $EXPID.*.data +/bin/rm $EXPID.*.ctl + +# Create DDF File +# --------------- +set monthlies = `/bin/ls -1 $EXPID.${TEM_Collection}.monthly.*.nc4` +set num = `/bin/ls -1 $monthlies | wc -l` +echo num = $num + + set edate = `echo $monthlies[$num] | cut -d "." -f${date_node}` + set eyear = `echo $edate | cut -c1-4` + set emonth = `echo $edate | cut -c5-6` + + set date = `echo $monthlies[1] | cut -d "." -f${date_node}` + set year = `echo $date | cut -c1-4` + set month = `echo $date | cut -c5-6` + +echo edate = $edate +echo eyear = $eyear +echo emonth = $emonth + +# set expdsc = `grep EXPDSC: $HISTORYRC | cut -d" " -f2` + set timinc = "mo" + + if( $month == "01" ) set MON = JAN + if( $month == "02" ) set MON = FEB + if( $month == "03" ) set MON = MAR + if( $month == "04" ) set MON = APR + if( $month == "05" ) set MON = MAY + if( $month == "06" ) set MON = JUN + if( $month == "07" ) set MON = JUL + if( $month == "08" ) set MON = AUG + if( $month == "09" ) set MON = SEP + if( $month == "10" ) set MON = OCT + if( $month == "11" ) set MON = NOV + if( $month == "12" ) set MON = DEC + +@ nmonths = ( $eyear - $year ) * 12 + ( $emonth - $month ) + 1 + + /bin/rm -f $EXPID.${TEM_Collection}.monthly.ddf + echo DSET ^$EXPID.${TEM_Collection}.monthly.%y4%m201.nc4 > $EXPID.${TEM_Collection}.monthly.ddf + echo TITLE Dummy Experiment Description >> $EXPID.${TEM_Collection}.monthly.ddf + echo OPTIONS template >> $EXPID.${TEM_Collection}.monthly.ddf + echo TDEF time $nmonths LINEAR 00:00Z01$MON$year 1$timinc >> $EXPID.${TEM_Collection}.monthly.ddf + + + + + diff --git a/plots/grads_util/getinfo b/plots/grads_util/getinfo index 5aa2a668..cf6a73ea 100644 --- a/plots/grads_util/getinfo +++ b/plots/grads_util/getinfo @@ -17,6 +17,19 @@ grav = 9.80665 * ----------------- +* Get Grads Config and Version +* ---------------------------- +if (name="config") + 'query config' + config = sublin(result,1) +endif + +if (name="version") + 'query config' + config = sublin(result,1) + version = subwrd(config,2) +endif + if( file = "" ) 'query file' file = sublin(result,1) @@ -318,6 +331,9 @@ pagey = subwrd(pagey ,6) * Return Options * -------------- + if (name="config" ) ; return config ; endif + if (name="version" ) ; return version ; endif + if (name="dset" ) ; return dset ; endif if (name="file" ) ; return file ; endif if (name="desc" ) ; return desc ; endif diff --git a/plots/grads_util/getobs.gs b/plots/grads_util/getobs.gs index 60caa711..922c1252 100755 --- a/plots/grads_util/getobs.gs +++ b/plots/grads_util/getobs.gs @@ -42,6 +42,7 @@ endif if( file = "NULL" ) if( ftype = "BIN" ) ; ' open 'fname ; endif if( ftype = "HDF" ) ; 'xdfopen 'fname ; endif +if( ftype = "CFIO" ) ; 'xdfopen 'fname ; endif if( ftype = "SDF" ) ; 'sdfopen 'fname ; endif diff --git a/plots/grads_util/makezdif.gs b/plots/grads_util/makezdif.gs index bd9abf0b..c8aad89c 100644 --- a/plots/grads_util/makezdif.gs +++ b/plots/grads_util/makezdif.gs @@ -359,9 +359,9 @@ minval = 1e15 endwhile *say ' ' -*'close 'newfile '!remove ZDIFILE.txt' -'run setenv "ZDIFILE" 'newfile +say 'run setenv "ZDIFILE" 'newfile + 'run setenv "ZDIFILE" 'newfile * Reset Initial Environment Settings * ---------------------------------- diff --git a/plots/grads_util/rgbset.gs b/plots/grads_util/rgbset.gs index 872d3f57..fc195420 100644 --- a/plots/grads_util/rgbset.gs +++ b/plots/grads_util/rgbset.gs @@ -90,6 +90,13 @@ 'set rgb 93 234 234 234' 'set rgb 94 252 252 252' +* RGB Parameters for STATS Package using Latest Grads Version +* ----------------------------------------------------------- +'getinfo version' + version = result + +if( version != 'v2.0.1.oga.1' ) + * black dots 'set tile 200 2 10 10 2 0' 'set rgb 200 tile 200' @@ -181,3 +188,5 @@ 'set tile 229 2 15 15 1 29' 'set rgb 229 tile 229' +endif + diff --git a/plots/grads_util/writegrads.gs b/plots/grads_util/writegrads.gs index d425b35e..a3822a35 100755 --- a/plots/grads_util/writegrads.gs +++ b/plots/grads_util/writegrads.gs @@ -1,11 +1,16 @@ function writegrads(args) ***************************************************************************************** -* Note: To Write Data in Reverse Order -* run setenv ZFLIP ON +* * args: -vars list of variables to write (Default: all variables) * -levs list of levels to write (Default: all levels) +* -time all (all times to single file, Default) +* one (one time per file) * -name name of binary name.data and name.ctl file (Default: name = grads) +* +* Note: To Write Data in Reverse Order +* run setenv ZFLIP ON +* ***************************************************************************************** 'numargs 'args @@ -36,6 +41,7 @@ if( xdim = 1 ) ; dlon = 1 ; endif VARS = 'ALL' LEVS = 'ALL' +TIME = 'ALL' NAME = 'grads' n = 0 @@ -43,6 +49,7 @@ NAME = 'grads' while ( num < numargs ) num = num + 1 if( subwrd(args,num) = '-name' ) ; NAME = subwrd(args,num+1) ; endif + if( subwrd(args,num) = '-time' ) ; TIME = subwrd(args,num+1) ; endif * Read VARS * --------- @@ -121,6 +128,9 @@ endif 'numargs 'VARS numvars = result +'run uppercase 'TIME + TIME = result + ******************************************* 'q gxout' @@ -139,14 +149,6 @@ endif lonbeg = result 'getinfo lat' latbeg = result -'setx' -'sety' - -* -------------------------------- - -'set t 1' -'getinfo date' - begdate = result 'getinfo tinc' tinc = result 'getinfo tunit' @@ -159,29 +161,53 @@ if( tunit = 'month' ) ; tunit = mo ; endif if( tunit = 'day' ) ; tunit = dy ; endif if( tunit = 'hour' ) ; tunit = hr ; endif +'setx' +'sety' + +* -------------------------------- + +* -------------------------------------------------------------------- + +t=1 +while(t<=tdim) + + 'set t 't + +'getinfo date' + date = result + year = substr( date,9,4 ) + day = substr( date,4,2 ) + month = get_month() + yyyymmdd = year''month''day + +if( t = 1 | TIME = 'ONE' ) + * Write 'NAME'.ctl * ---------------- -'!remove 'NAME'.ctl' -'!touch 'NAME'.ctl' -'!echo dset ^'NAME'.data >> 'NAME'.ctl' -'!echo title 'label' >> 'NAME'.ctl' -'!echo undef 1e15 >> 'NAME'.ctl' -'!echo xdef 'xdim' linear 'lonbeg' 'dlon' >> 'NAME'.ctl' -'!echo ydef 'ydim' linear 'latbeg' 'dlat' >> 'NAME'.ctl' -'!echo zdef 'zdim' levels 'LEVS' >> 'NAME'.ctl' -'!echo tdef 'tdim' linear 'begdate' 'tinc''tunit' >> 'NAME'.ctl' -'!echo vars 'numvars' >> 'NAME'.ctl' +'!remove 'NAME'.'yyyymmdd'.ctl' +'!touch 'NAME'.'yyyymmdd'.ctl' +'!echo dset ^'NAME'.'yyyymmdd'.data >> 'NAME'.'yyyymmdd'.ctl' +'!echo title 'label' >> 'NAME'.'yyyymmdd'.ctl' +'!echo undef 1e15 >> 'NAME'.'yyyymmdd'.ctl' +'!echo xdef 'xdim' linear 'lonbeg' 'dlon' >> 'NAME'.'yyyymmdd'.ctl' +'!echo ydef 'ydim' linear 'latbeg' 'dlat' >> 'NAME'.'yyyymmdd'.ctl' +'!echo zdef 'zdim' levels 'LEVS' >> 'NAME'.'yyyymmdd'.ctl' +if( TIME = 'ALL' ) +'!echo tdef 'tdim' linear 'date' 'tinc''tunit' >> 'NAME'.'yyyymmdd'.ctl' +else +'!echo tdef 1 linear 'date' 'tinc''tunit' >> 'NAME'.'yyyymmdd'.ctl' +endif +'!echo vars 'numvars' >> 'NAME'.'yyyymmdd'.ctl' * -------------------------------------------------------------------- 'set gxout fwrite' -'set fwrite 'NAME'.data' +'set fwrite 'NAME'.'yyyymmdd'.data' -'set undef 1e15' +endif -t=1 -while(t<=tdim) - 'set t 't + +'set undef 1e15' n=1 while(n<=nvars) @@ -232,7 +258,7 @@ while(n<=nvars) if( name = slp ) ; zref = 0 ; endif if( zref = 0 ) - if( t=1 ) ; '!echo "'desc'" >> 'NAME'.ctl' ; endif + if( ( t=1 & TIME = 'ALL' ) | TIME = 'ONE' ) ; '!echo "'desc'" >> 'NAME'.'yyyymmdd'.ctl' ; endif * say 'Writing Variable: 'name 'set z 1' e = 1 @@ -242,7 +268,7 @@ while(n<=nvars) e = e + 1 endwhile else - if( t=1 ) ; '!echo "'desc'" >> 'NAME'.ctl' ; endif + if( ( t=1 & TIME = 'ALL' ) | TIME = 'ONE' ) ; '!echo "'desc'" >> 'NAME'.'yyyymmdd'.ctl' ; endif if(zflip != 'ON' ) z=1 @@ -286,11 +312,32 @@ while(n<=nvars) n=n+1 endwhile -if( t=1 ) ; '!echo endvars >> 'NAME'.ctl' ; endif + +if( ( t=1 & TIME = 'ALL' ) | TIME = 'ONE' ) ; '!echo endvars >> 'NAME'.'yyyymmdd'.ctl' ; endif + if( TIME = 'ONE' ) ; 'disable fwrite' ; 'set gxout 'gxout ; endif + say ' ' t=t+1 endwhile -'disable fwrite' +if( TIME = 'ALL' ) ; 'disable fwrite' ; endif 'set gxout 'gxout return + +function get_month(args) +'getinfo month' + month = result + if(month="JAN") ; monthnum = 01 ; endif + if(month="FEB") ; monthnum = 02 ; endif + if(month="MAR") ; monthnum = 03 ; endif + if(month="APR") ; monthnum = 04 ; endif + if(month="MAY") ; monthnum = 05 ; endif + if(month="JUN") ; monthnum = 06 ; endif + if(month="JUL") ; monthnum = 07 ; endif + if(month="AUG") ; monthnum = 08 ; endif + if(month="SEP") ; monthnum = 09 ; endif + if(month="OCT") ; monthnum = 10 ; endif + if(month="NOV") ; monthnum = 11 ; endif + if(month="DEC") ; monthnum = 12 ; endif +return monthnum + diff --git a/plots/hcmp/hcmp2 b/plots/hcmp/hcmp2 index 3a1c4135..17b16a7b 100644 --- a/plots/hcmp/hcmp2 +++ b/plots/hcmp/hcmp2 @@ -691,6 +691,8 @@ endwhile 'define qcmp'numexp''season' = sqrt(qcmp'numexp''season')' endif + 'define zmod'season' = regrid2( qmod'season' ,0.625,0.5,bs_p1,0,-90 )' + 'define zcmp'numexp''season' = regrid2( qcmp'numexp''season',0.625,0.5,bs_p1,0,-90 )' 'run setenv "CLIMATE" 'climcmp.ananam flag = "" @@ -756,7 +758,9 @@ say 'Looping through Comparison Experiments for Closeness Plots' say '----------------------------------------------------------' k = 1 while( k <= numexp ) -say 'Checking Comparison Experiment k = 'k' CTAG = 'ctag.k' TYPE = 'type.k' for Closeness to Verifications (MERRA-2 ... CMPEXP:V)' + say 'k = 'k' numexp = 'numexp + say 'Checking Comparison Experiment k = 'k' CTAG = 'ctag.k' TYPE = 'type.k' for Closeness to Verifications (MERRA-2 ... CMPEXP:V)' + pause TAG = k @@ -785,7 +789,9 @@ say 'Checking Comparison Experiment k = 'k' CTAG = 'ctag.k' TYPE = 'type.k' for * -------------------------------------------------------------------- n = 1 while( n <= numexp ) - say 'Testing 'qtag' and 'ctag.n' for closeness with 'ctag.TAG + say 'n = 'n' numexp = 'numexp + say 'Testing 'qtag' and 'ctag.n' for closeness with 'ctag.TAG + pause if( ctag.n != "merra" & ctag.n != "MERRA-2" & ctag.n != ctag.TAG & type.n != V & cname.n != 'NULL' ) say 'Closeness plot between exp: 'qtag @@ -800,12 +806,11 @@ say 'Checking Comparison Experiment k = 'k' CTAG = 'ctag.k' TYPE = 'type.k' for say '' 'define zcmp'TAG''season' = regrid2( qcmp'TAG''season',0.625,0.5,bs_p1,0,-90 )' - 'define zcmp'n''season' = regrid2( qcmp'n''season' ,0.625,0.5,bs_p1,0,-90 )' - 'define zmod'season' = regrid2( qmod'season' ,0.625,0.5,bs_p1,0,-90 )' flag = "" while ( flag = "" ) say 'Performing Closeness plots to: 'ctag.TAG' k = 'k + pause 'closeness -CVAR 'zcmp''n' -MVAR 'zmod' -OVAR 'zcmp''TAG' -CNAME 'ctag.n' -MNAME 'EXPORT' -ONAME 'ctag.TAG' -CDESC 'cdesc.n' -MDESC 'qdesc' -ODESC 'cdesc.TAG' -MFILE 'qfile' -MBEGDATE 'bdate' -MENDDATE 'edate' -CFILE 'cfile.n' -OFILE 'cfile.TAG' -CBEGDATE 'begdate.n' -CENDDATE 'enddate.n' -OBEGDATE 'begdate.k' -OENDDATE 'enddate.k' -EXPID 'expid' -PREFIX 'NULL' -SEASON 'season' -OUTPUT 'output' -CLIMEXP 'climexp' -CLIMCMP 'climcmp.nname' -CLIMOBS 'climcmp.kname' -GC 'GC' -MATH 'NULL' -LEVEL 'level * 'closeness -CVAR 'zcmp''n @@ -863,6 +868,7 @@ say 'Checking Comparison Experiment k = 'k' CTAG = 'ctag.k' TYPE = 'type.k' for say ' ' say 'Closeness of 'cname.k' of type: 'type.k' to 'ctag.k' NOT Performed' say '------------------------------------------------------------------' + pause say ' ' endif ;* END ctag.k Closeness If Test @@ -1332,6 +1338,7 @@ say ' climcmp: 'climcmp.cname say ' climobs: 'climobs.ananam say ' Total numexp: 'numexp say '' +pause flag = "" while ( flag = "" ) diff --git a/plots/portrait.script b/plots/portrait.script index 24eb7f16..b32b1834 100755 --- a/plots/portrait.script +++ b/plots/portrait.script @@ -279,7 +279,8 @@ endif # Residual Circulation -# -------------------- +# Note: Always run MAKERES before running the Residual PORTRAIT and LANDSCAPE plots +# --------------------------------------------------------------------------------- if( $plot == residual ) then $GEOSUTIL/plots/res/makeres $source @@ -289,7 +290,7 @@ if( $plot == residual ) then setenv ORIENTATION LANDSCAPE $grads $batch -l -c "run $GEOSUTIL/plots/res/res $expid $output $debug $seasons " $grads $batch -l -c "run $GEOSUTIL/plots/res/setup_wstar $source $expid $output $seasons " - $grads $batch -l -c "run $GEOSUTIL/plots/res/setup_epflx $source $expid $output $seasons " + $grads $batch -l -c "run $GEOSUTIL/plots/res/setup_epflx $source $expid $output $debug $seasons " $grads $batch -l -c "run $GEOSUTIL/plots/zcmp/zcmp $expid $output $seasons " endif diff --git a/plots/quickplot b/plots/quickplot index 18f79615..b7ab378c 100755 --- a/plots/quickplot +++ b/plots/quickplot @@ -353,14 +353,12 @@ source $plots_dir/$plot.$datetime.$type/.quickplotrc # Set Default OPS Comparison Experiments # -------------------------------------- -set cmpops = "$VERIFICATION/MERRA2_MEANS:A \ - $VERIFICATION/ERA5_Monthly:V " +# set cmpops = "$VERIFICATION/MERRA2_MEANS:A " +# set cmpops = "$VERIFICATION/ERA5_Monthly:V " + set cmpops = "$VERIFICATION/MERRA2_MEANS:A $VERIFICATION/ERA5_Monthly:V " +# set cmpops = " " -# $VERIFICATION/MERRA_MEANS:A \ -#set cmpops = "$VERIFICATION/ERA5_Monthly:V " -#set cmpops = " " - - set cmpops = `echo $cmpops` + set cmpops = `echo $cmpops` # Update CMPEXP # ------------- @@ -393,6 +391,7 @@ echo " ANALYSIS = " $ANALYSIS echo " STD_DEV = " $STD_DEV echo " PLOT = " $plot echo " CMPEXP = " $CMPEXP +echo " CMPOPS = " $cmpops echo " CMPEXP_ONLY = " $cmpexp_only echo " CMPOPS_ONLY = " $cmpops_only echo " " diff --git a/plots/res/VERIFICATION.MERRA-2.rc b/plots/res/VERIFICATION.MERRA-2.rc index b28a8fa4..e1ace0c5 100644 --- a/plots/res/VERIFICATION.MERRA-2.rc +++ b/plots/res/VERIFICATION.MERRA-2.rc @@ -1,8 +1,8 @@ OBSNAM: MERRA-2 OBSDSC: MERRA-2_Reanalysis - filename: 'VERIFICATION/MERRA2_MEANS/residual.MERRA-2.ctl' - filetype: 'BIN' + filename: 'VERIFICATION/MERRA2_MEANS/inst3_3d_tem_Np/MERRA-2.inst3_3d_tem_Np.monthly.ddf' + filetype: 'CFIO' fields: 'STR' , 'DYN' , 'STR' , 1 'RES' , 'DYN' , 'RES' , 1 'DUM' , 'DYN' , 'DUM' , 1 diff --git a/plots/res/VERIFICATION.rc.tmpl b/plots/res/VERIFICATION.rc.tmpl index 41c92deb..2bd43fb1 100644 --- a/plots/res/VERIFICATION.rc.tmpl +++ b/plots/res/VERIFICATION.rc.tmpl @@ -1,8 +1,8 @@ OBSNAM: @EXPID OBSDSC: @EXPDSC - filename: '@pwd/residual.@EXPID.ctl' - filetype: 'BIN' + filename: '@filename' + filetype: 'CFIO' fields: 'STR' , 'DYN' , 'STR' , 1 'RES' , 'DYN' , 'RES' , 1 'DUM' , 'DYN' , 'DUM' , 1 diff --git a/plots/res/epflx.gs b/plots/res/epflx.gs index e1f7e682..b1346463 100644 --- a/plots/res/epflx.gs +++ b/plots/res/epflx.gs @@ -1,19 +1,32 @@ function epflx (args) -expid = subwrd(args,1) -season = subwrd(args,2) -index = subwrd(args,3) -output = subwrd(args,4) + expid = subwrd(args,1) +TEM_Collection = subwrd(args,2) + season = subwrd(args,3) + index = subwrd(args,4) + output = subwrd(args,5) time = substr(index,1,1) num = substr(index,2,2) +say 'Inside epflx, expid = 'expid +say 'Inside epflx, season = 'season +say 'Inside epflx, index = 'index +say 'Inside epflx, output = 'output +say 'Inside epflx, time = 'time +say 'Inside epflx, dfile = 'num +pause 'set dfile 'num -'setz ' +'q file' +say result +pause + + 'setz ' 'getinfo desc' desc = result + say 'desc = 'desc node = '' - m = 2 + m = 1 '!remove NODE.txt' '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' @@ -25,7 +38,11 @@ num = substr(index,2,2) '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' node = result - while( node != 'ctl' ) + + say 'NODE = 'node + say 'TEM_Collection = 'TEM_Collection + + while( node != TEM_Collection ) EXP = EXP'.'node say 'EXP = 'EXP m = m + 1 @@ -35,15 +52,19 @@ num = substr(index,2,2) node = result endwhile desc = EXP + say 'EXP: 'EXP + pause 'run getenv MASKFILE' maskfile = result if( time = 'A' ) + 'setdates' 'run getenv BEGDATE' begdate = result 'run getenv ENDDATE' enddate = result + 'set t 1' else 'set t 1' 'run getinfo date' @@ -86,6 +107,9 @@ say 'ENDDATE = 'enddate ' set arrlab off ' ' set ylopts 1 3 0.15 ' +* ------------------------------------------------------------------------------------- +* ------------------------------------------------------------------------------------- + ' vpage 1 1 2 2 -top 0.6 -bot -0.10' ' set ccolor rainbow' ' set csmooth on' @@ -111,6 +135,8 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' set strsiz 0.10' ' draw string 8.25 0.7 Z/Y Ratio: 'arrfct +* ------------------------------------------------------------------------------------- +* ------------------------------------------------------------------------------------- ' vpage 2 1 2 2 -top 0.6 -bot -0.10' ' set ccolor rainbow' @@ -136,6 +162,9 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' set strsiz 0.10' ' draw string 8.25 0.7 Ratio Z/Y: 'arrfct +* ------------------------------------------------------------------------------------- +* ------------------------------------------------------------------------------------- + ' vpage 1 2 2 2 -top 0.25 -bot 0.25' ' set ccolor rainbow' ' set csmooth on' @@ -160,6 +189,9 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' set strsiz 0.10' ' draw string 8.25 0.7 Ratio Z/Y: 'arrfct +* ------------------------------------------------------------------------------------- +* ------------------------------------------------------------------------------------- + ' vpage 2 2 2 2 -top 0.25 -bot 0.25' ' set csmooth on' ' set ccolor rainbow' @@ -186,6 +218,8 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' draw string 8.25 0.7 Ratio Z/Y: 'arrfct +* ------------------------------------------------------------------------------------- +* ------------------------------------------------------------------------------------- ' set vpage off ' diff --git a/plots/res/epflx_diff.gs b/plots/res/epflx_diff.gs index d4fa3b91..63e4b23e 100644 --- a/plots/res/epflx_diff.gs +++ b/plots/res/epflx_diff.gs @@ -2,17 +2,22 @@ function epflx_diff (args) expid = subwrd(args,1) cmpid = subwrd(args,2) -season = subwrd(args,3) -index = subwrd(args,4) -index2 = subwrd(args,5) -output = subwrd(args,6) -detail = subwrd(args,7) - +TEM_expid = subwrd(args,3) +TEM_cmpid = subwrd(args,4) +season = subwrd(args,5) +index = subwrd(args,6) +index2 = subwrd(args,7) +output = subwrd(args,8) +detail = subwrd(args,9) +detail = d time = substr(index ,1,1) num = substr(index ,2,2) num2 = substr(index2,2,2) +*'run getenv "TEM_Collection"' +* TEM_Collection = result +*say 'TEM_Collection = 'TEM_Collection 'set dfile 'num2 'setz ' @@ -20,7 +25,7 @@ num2 = substr(index2,2,2) desc = result node = '' - m = 2 + m = 1 '!remove NODE.txt' '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' @@ -32,7 +37,7 @@ num2 = substr(index2,2,2) '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' node = result - while( node != 'ctl' ) + while( node != TEM_cmpid ) EXP = EXP'.'node say 'EXP = 'EXP m = m + 1 @@ -50,7 +55,7 @@ num2 = substr(index2,2,2) desc = result node = '' - m = 2 + m = 1 '!remove NODE.txt' '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' @@ -62,7 +67,7 @@ num2 = substr(index2,2,2) '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' node = result - while( node != 'ctl' ) + while( node != TEM_expid ) EXP = EXP'.'node say 'EXP = 'EXP m = m + 1 @@ -77,10 +82,12 @@ num2 = substr(index2,2,2) maskfile = result if( time = 'A' ) + 'setdates' 'run getenv BEGDATE' begdate = result 'run getenv ENDDATE' enddate = result + 'set t 1' else 'set t 1' 'run getinfo date' @@ -134,6 +141,9 @@ say 'ENDDATE = 'enddate 'set grads off ' 'set ylopts 1 3 0.15 ' +* ---------------------------------------------------------------------- +* ---------------------------------------------------------------------- + ' vpage 1 1 2 2 -top 0.6 -bot -0.10' ' set lev 1000 0.8' ' set lev 925 0.8' @@ -143,30 +153,71 @@ say 'ENDDATE = 'enddate ' d epdiff' ' cbarn -xmid 6 -ndot 0' ' draw ylab hPa ' -' set ccolor 1' -if( detail != '' ) - arrlen = 1 - arrscl = 1e7 - arrfct = 5 -else - arrlen = 1 - arrscl = 1e8 - arrfct = 50 -endif - 'define arry = epydiffz' - 'define arrz = -epzdiffz*'arrfct + +*if( detail != '' ) +* arrlen = 1 +* arrscl = 1e7 +* arrfct = 5 +* arrscl = 5e6 +* arrfct = 10 +*else +* arrlen = 1 +* arrscl = 1e8 +* arrfct = 50 +*endif + + 'qminmax epydiffz/1e6' + epydiffmin = subwrd(result,1) + epydiffmax = subwrd(result,2) + epydiffmag = epydiffmax - epydiffmin + say 'epydiffmax = 'epydiffmax + say 'epydiffmin = 'epydiffmin + say 'epydiffmag = 'epydiffmag + say ' ' + + 'qminmax epzdiffz/1e6' + epzdiffmin = subwrd(result,1) + epzdiffmax = subwrd(result,2) + epzdiffmag = epzdiffmax - epzdiffmin + say 'epzdiffmax = 'epzdiffmax + say 'epzdiffmin = 'epzdiffmin + say 'epzdiffmag = 'epzdiffmag + say ' ' + + epyzratio = epydiffmag / epzdiffmag + epyzratio = epyzratio * 100 + 'getint 'epyzratio + arrfct = result / 100 + say 'epyzratio = 'arrfct + say ' ' + + 'define arry = epydiffz/1e6' + 'define arrz = -epzdiffz/1e6*'arrfct + 'define arrmag = sqrt( arry*arry + arrz*arrz )' + + 'qminmax arrmag' + arrmag = subwrd(result,2) + 'define arryskp = skip(arry,'skipval')' + arrlen = 1 + arrscl = arrmag 'set arrscl 'arrlen' 'arrscl + 'set ccolor 1' 'd arryskp;arrz' -xrit = 4.0 -ybot = 0.6 -rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + + xrit = 4.0 + ybot = 0.6 + rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + ' set string 1 c 4 ' ' set strsiz 0.09' ' draw string 8.25 0.7 Z/Y Ratio: 'arrfct +* ---------------------------------------------------------------------- +* ---------------------------------------------------------------------- ' vpage 2 1 2 2 -top 0.6 -bot -0.10' +' set gxout shaded' ' set csmooth on' ' set lev 1000 8' ' set lev 925 8' @@ -176,29 +227,71 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' cbarn -xmid 6 -ndot 0' ' draw ylab hPa ' ' set ccolor 1' -if( detail != '' ) - arrlen = 1 - arrscl = 1e7 - arrfct = 5 -else - arrlen = 1 - arrscl = 1e8 - arrfct = 50 -endif - 'define arry = epydiffz' - 'define arrz = -epzdiffz*'arrfct + +*if( detail != '' ) +* arrlen = 1 +* arrscl = 1e7 +* arrfct = 5 +* arrscl = 5e6 +* arrfct = 10 +*else +* arrlen = 1 +* arrscl = 1e8 +* arrfct = 50 +*endif + + 'qminmax epydiffz/1e6' + epydiffmin = subwrd(result,1) + epydiffmax = subwrd(result,2) + epydiffmag = epydiffmax - epydiffmin + say 'epydiffmax = 'epydiffmax + say 'epydiffmin = 'epydiffmin + say 'epydiffmag = 'epydiffmag + say ' ' + + 'qminmax epzdiffz/1e6' + epzdiffmin = subwrd(result,1) + epzdiffmax = subwrd(result,2) + epzdiffmag = epzdiffmax - epzdiffmin + say 'epzdiffmax = 'epzdiffmax + say 'epzdiffmin = 'epzdiffmin + say 'epzdiffmag = 'epzdiffmag + say ' ' + + epyzratio = epydiffmag / epzdiffmag + epyzratio = epyzratio * 100 + 'getint 'epyzratio + arrfct = result / 100 + say 'epyzratio = 'arrfct + say ' ' + + 'define arry = epydiffz/1e6' + 'define arrz = -epzdiffz/1e6*'arrfct + 'define arrmag = sqrt( arry*arry + arrz*arrz )' + + 'qminmax arrmag' + arrmag = subwrd(result,2) + 'define arryskp = skip(arry,'skipval')' + arrlen = 1 + arrscl = arrmag 'set arrscl 'arrlen' 'arrscl + 'set ccolor 1' 'd arryskp;arrz' -xrit = 4.0 -ybot = 0.6 -rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + + xrit = 4.0 + ybot = 0.6 + rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + ' set string 1 c 4 ' ' set strsiz 0.09' ' draw string 8.25 0.7 Z/Y Ratio: 'arrfct +* ---------------------------------------------------------------------- +* ---------------------------------------------------------------------- ' vpage 1 2 2 2 -top 0.25 -bot 0.25' +' set gxout shaded' ' set csmooth on' ' set lev 1000 80' ' set lev 925 80' @@ -208,29 +301,71 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' cbarn -xmid 6 -ndot 0' ' draw ylab hPa ' ' set ccolor 1' -if( detail != '' ) - arrlen = 1 - arrscl = 1e7 - arrfct = 5 -else - arrlen = 1 - arrscl = 1e8 - arrfct = 50 -endif - 'define arry = epydiffz' - 'define arrz = -epzdiffz*'arrfct + +*if( detail != '' ) +* arrlen = 1 +* arrscl = 1e7 +* arrfct = 5 +* arrscl = 5e6 +* arrfct = 10 +*else +* arrlen = 1 +* arrscl = 1e8 +* arrfct = 50 +*endif + + 'qminmax epydiffz/1e6' + epydiffmin = subwrd(result,1) + epydiffmax = subwrd(result,2) + epydiffmag = epydiffmax - epydiffmin + say 'epydiffmax = 'epydiffmax + say 'epydiffmin = 'epydiffmin + say 'epydiffmag = 'epydiffmag + say ' ' + + 'qminmax epzdiffz/1e6' + epzdiffmin = subwrd(result,1) + epzdiffmax = subwrd(result,2) + epzdiffmag = epzdiffmax - epzdiffmin + say 'epzdiffmax = 'epzdiffmax + say 'epzdiffmin = 'epzdiffmin + say 'epzdiffmag = 'epzdiffmag + say ' ' + + epyzratio = epydiffmag / epzdiffmag + epyzratio = epyzratio * 100 + 'getint 'epyzratio + arrfct = result / 100 + say 'epyzratio = 'arrfct + say ' ' + + 'define arry = epydiffz/1e6' + 'define arrz = -epzdiffz/1e6*'arrfct + 'define arrmag = sqrt( arry*arry + arrz*arrz )' + + 'qminmax arrmag' + arrmag = subwrd(result,2) + 'define arryskp = skip(arry,'skipval')' + arrlen = 1 + arrscl = arrmag 'set arrscl 'arrlen' 'arrscl + 'set ccolor 1' 'd arryskp;arrz' -xrit = 4.0 -ybot = 0.6 -rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + + xrit = 4.0 + ybot = 0.6 + rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + ' set string 1 c 4 ' ' set strsiz 0.09' ' draw string 8.25 0.7 Z/Y Ratio: 'arrfct +* ---------------------------------------------------------------------- +* ---------------------------------------------------------------------- ' vpage 2 2 2 2 -top 0.25 -bot 0.25' +' set gxout shaded' ' set csmooth on' ' set grads off ' ' set lev 100 18' @@ -249,30 +384,68 @@ rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) ' cbarn -xmid 6 -ndot 1' ' draw ylab hPa ' ' set ccolor 1' -if( detail != '' ) - arrlen = 1 - arrscl = 5e6 - arrfct = 500 -else - arrlen = 1 - arrscl = 1e8 - arrfct = 500 -endif - - 'define arry = epydiffz' - 'define arrz = -epzdiffz*'arrfct +*if( detail != '' ) +* arrlen = 1 +* arrscl = 5e6 +* arrfct = 500 +* arrscl = 5e5 +* arrfct = 300 +*else +* arrlen = 1 +* arrscl = 1e8 +* arrfct = 500 +*endif + + 'qminmax epydiffz/1e6' + epydiffmin = subwrd(result,1) + epydiffmax = subwrd(result,2) + epydiffmag = epydiffmax - epydiffmin + say 'epydiffmax = 'epydiffmax + say 'epydiffmin = 'epydiffmin + say 'epydiffmag = 'epydiffmag + say ' ' + + 'qminmax epzdiffz/1e6' + epzdiffmin = subwrd(result,1) + epzdiffmax = subwrd(result,2) + epzdiffmag = epzdiffmax - epzdiffmin + say 'epzdiffmax = 'epzdiffmax + say 'epzdiffmin = 'epzdiffmin + say 'epzdiffmag = 'epzdiffmag + say ' ' + + epyzratio = epydiffmag / epzdiffmag + epyzratio = epyzratio * 100 + 'getint 'epyzratio + arrfct = result / 100 + say 'epyzratio = 'arrfct + say ' ' + + 'define arry = epydiffz/1e6' + 'define arrz = -epzdiffz/1e6*'arrfct + 'define arrmag = sqrt( arry*arry + arrz*arrz )' + + 'qminmax arrmag' + arrmag = subwrd(result,2) + 'define arryskp = skip(arry,'skipval')' + arrlen = 1 + arrscl = arrmag 'set arrscl 'arrlen' 'arrscl + 'set ccolor 1' 'd arryskp;arrz' -xrit = 4.0 -ybot = 0.6 -rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + + xrit = 4.0 + ybot = 0.6 + rc = arrow(xrit-0.25,ybot+0.2,arrlen,arrscl) + ' set string 1 c 4 ' ' set strsiz 0.09' ' draw string 8.25 0.7 Z/Y Ratio: 'arrfct - +* ---------------------------------------------------------------------- +* ---------------------------------------------------------------------- ' set vpage off ' diff --git a/plots/res/plot_season.gs b/plots/res/plot_season.gs index fad0ea39..5e12daa8 100644 --- a/plots/res/plot_season.gs +++ b/plots/res/plot_season.gs @@ -129,13 +129,56 @@ endif * ------------------------------------------------------------------------ * ------------------------------------------------------------------------ +* Loop over Experiment Datasets to get Experiment IDs +* --------------------------------------------------- +say '' +'getpwd' + pwd = result + +'getnumrc 'pwd + rcinfo = result + numrc = subwrd( rcinfo,1 ) + num = 1 + cnt = 0 +while( num <= numrc ) + loc = num + 1 + rcfile = subwrd( rcinfo,loc ) + say 'num = 'num' rcfile = 'rcfile + + '!grep filename 'rcfile' | cut -d'"\' -f2 > CTLFILE.txt" + 'run getenv CTLFILE ' + CTLFILE.num = result + + '!basename 'rcfile' > BASENAME.txt' + 'run getenv BASENAME' + basename = result + length = strlen(basename) - 3 + say 'original length = 'length + '!echo 'basename' | cut -b1-'length' | cut -d. -f2 > CMPID.txt' + 'run getenv CMPID ' + CMPID.num = result + + say ' CMPID #'num' = 'CMPID.num + say ' CTLFILE #'num' = 'CTLFILE.num + + '!remove TEM_NAME.txt.' + '!basename 'CTLFILE.num' | cut -d. -f2 > TEM_NAME.txt' + say 'run getenv "TEM_NAME"' + 'run getenv "TEM_NAME"' + TEM_NAME.num = result + say 'TEM_Collection = 'TEM_NAME.num +* pause + + num = num + 1 +endwhile + n = 1 while( n<=numfiles ) 'set dfile 'n 'getinfo desc' desc = result node = '' - m = 2 + m = 1 '!remove NODE.txt' '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' @@ -147,7 +190,13 @@ endif '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' node = result - while( node != 'ctl' & node != 'data' ) + + TEM_Collection = TEM_NAME.n + say 'NODE = 'node + say 'TEM_Collection = 'TEM_Collection + + while( node != TEM_Collection ) + EXP.n = EXP.n'.'node say 'EXP'n' = 'EXP.n m = m + 1 @@ -477,8 +526,8 @@ endwhile say 'Final axmin: 'axmin' axmax: 'axmax say ' ' - axmin = 0.21509 - axmax = 0.44501 +* axmin = 0.21509 +* axmax = 0.44501 'set axlim 'axmin' 'axmax id.1 = 'a' @@ -926,10 +975,10 @@ endif 'set parea off' if( tags = '' ) - filename = 'WSTAR_Profile_using_'talatsz'_TALATS.'season + filename = 'WSTAR_using_'talatsz'_TALATS.'season else -* filename = 'WSTAR_Profile_using_'talatsz'_TALATS.'tags'.'season - filename = 'WSTAR_Profile_using_'talatsz'_TALATS.'exps'.'season +* filename = 'WSTAR_using_'talatsz'_TALATS.'tags'.'season + filename = 'WSTAR_using_'talatsz'_TALATS.'exps'.'season endif 'myprint -name 'output'/'filename diff --git a/plots/res/res b/plots/res/res index 9ae1f9d6..db77f517 100644 --- a/plots/res/res +++ b/plots/res/res @@ -56,6 +56,7 @@ endwhile 'set dfile 'modfil2 'sety' 'setz' +'getdates' 'seasonal res mod' @@ -88,11 +89,13 @@ if( obsid != expid ) 'set dfile 'obsfil1 'sety' 'setz' +'getdates' 'seasonal str obs' 'set dfile 'obsfil2 'sety' 'setz' +'getdates' 'seasonal res obs' 'sett' @@ -121,6 +124,7 @@ else 'count "'season'" 'begdate' 'enddate nmod = result + * Make LANDSCAPE Plots * -------------------- if( orientation = LANDSCAPE ) @@ -151,6 +155,9 @@ else flag = "" while ( flag = "" ) + say ' ' + say 'run 'geosutil'/plots/res/res_2 'expid' 'season' 'output' 'obsid' 'obsdsc' 'nmod' 'nobs' 'expdsc' 'begdateo' 'enddateo + say ' ' 'run 'geosutil'/plots/res/res_2 'expid' 'season' 'output' 'obsid' 'obsdsc' 'nmod' 'nobs' 'expdsc' 'begdateo' 'enddateo if( debug = "debug" ) say "Hit ENTER to repeat plot" @@ -164,6 +171,9 @@ else flag = "" while ( flag = "" ) + say ' ' + say 'run 'geosutil'/plots/res/res_3 'expid' 'season' 'output' 'obsid' 'obsdsc' 'nmod' 'nobs' 'expdsc' 'begdateo' 'enddateo + say ' ' 'run 'geosutil'/plots/res/res_3 'expid' 'season' 'output' 'obsid' 'obsdsc' 'nmod' 'nobs' 'expdsc' 'begdateo' 'enddateo if( debug = "debug" ) say "Hit ENTER to repeat plot" diff --git a/plots/res/res_2.gs b/plots/res/res_2.gs index 46fcb5a3..5da0676b 100644 --- a/plots/res/res_2.gs +++ b/plots/res/res_2.gs @@ -11,6 +11,9 @@ function res_2 (args) begdateo = subwrd(args,9) enddateo = subwrd(args,10) +'getinfo file' + dfile = result + * Get Dates for Plots * ------------------- 'run getenv "BEGDATE"' @@ -21,6 +24,9 @@ function res_2 (args) 'run getenv "CINTDIFF"' CINTDIFF = result +'run getenv "ZDIFILE" ' + zdifile = result + 'run getenv "CLIMATE"' climate = result if(begdate = begdateo & enddate = enddateo ) @@ -43,12 +49,6 @@ endif ******************************************************** -if( CINTDIFF != NULL ) -* -------------------- - -'getinfo file' - dfile = result - * Get Bounding Values * ------------------- 'getinfo xdim' @@ -106,19 +106,19 @@ z1 = result z2 = result endif -'run getenv "ZDIFILE" ' - zdifile = result -'set dfile ' zdifile -'q file' -say 'ZDIFILE: 'result -'set x 1' -'set t 1' -'sety' -'setz' - -'minmax strdifz' - dqmax = subwrd(result,1) - dqmin = subwrd(result,2) +if( CINTDIFF != NULL ) +* -------------------- + 'run getenv "ZDIFILE" ' + zdifile = result + 'set dfile ' zdifile + 'set x 1' + 'set t 1' + 'sety' + 'setz' + 'minmax strdifz' + dqmax = subwrd(result,1) + dqmin = subwrd(result,2) +endif * Reset initial space-time environment * ------------------------------------ @@ -128,8 +128,6 @@ say 'ZDIFILE: 'result 'set z 'z1' 'z2 'set t 1' -endif - ******************************************************** * Streamfunction ******************************************************** @@ -192,11 +190,13 @@ endif 'set t 1' 'd str'season'obs' -'set dfile ' zdifile +'set dfile 'zdifile 'set x 1' 'set t 1' 'sety' 'setz' +'getinfo undef' + undef = result if( CINTDIFF != NULL ) * -------------------- @@ -267,8 +267,19 @@ if( CINTDIFF != NULL ) else dn = 0 - 'shades 0.3' - 'd maskout( strdifz,abs(strdifz)-0.3 )' + 'define test = maskout( strdifz,abs(strdifz)-0.3 )' + 'minmax test' + testmax = subwrd(result,1) + testmin = subwrd(result,2) + if( testmax = undef & testmin = undef ) + 'shades strdifz 0' + cint = result + 'shades 'cint + 'd maskout( strdifz,abs(strdifz)-'cint' )' + else + 'shades 0.3' + 'd maskout( strdifz,abs(strdifz)-0.3 )' + endif endif @@ -287,7 +298,6 @@ endif 'set cmin -40' 'd str'season'obs' - 'set vpage off' 'set parea off' 'set vpage 0 8.5 0.0 11' diff --git a/plots/res/res_3.gs b/plots/res/res_3.gs index cbbe220f..96d5df30 100644 --- a/plots/res/res_3.gs +++ b/plots/res/res_3.gs @@ -1,4 +1,4 @@ -function res_2 (args) +function res_3 (args) expid = subwrd(args,1) season = subwrd(args,2) @@ -11,6 +11,9 @@ function res_2 (args) begdateo = subwrd(args,9) enddateo = subwrd(args,10) +'getinfo file' + dfile = result + * Get Dates for Plots * ------------------- 'run getenv "BEGDATE"' @@ -21,6 +24,9 @@ function res_2 (args) 'run getenv "CINTDIFF"' CINTDIFF = result +'run getenv "ZDIFILE" ' + zdifile = result + 'run getenv "CLIMATE"' climate = result if(begdate = begdateo & enddate = enddateo ) @@ -43,12 +49,6 @@ endif ******************************************************** -if( CINTDIFF != NULL ) -* -------------------- - -'getinfo file' - dfile = result - * Get Bounding Values * ------------------- 'getinfo xdim' @@ -106,20 +106,20 @@ z1 = result z2 = result endif -'run getenv "ZDIFILE" ' - zdifile = result -'set dfile ' zdifile -'q file' -say 'ZDIFILE: 'result -'set x 1' -'set t 1' -'sety' -'setz' - -'set lev 900 1' -'minmax resdifz' - dqmax = subwrd(result,1) - dqmin = subwrd(result,2) +if( CINTDIFF != NULL ) +* -------------------- + 'run getenv "ZDIFILE" ' + zdifile = result + 'set dfile ' zdifile + 'set x 1' + 'set t 1' + 'sety' + 'setz' + 'set lev 900 1' + 'minmax resdifz' + dqmax = subwrd(result,1) + dqmin = subwrd(result,2) +endif * Reset initial space-time environment * ------------------------------------ @@ -129,8 +129,6 @@ say 'ZDIFILE: 'result 'set z 'z1' 'z2 'set t 1' -endif - ******************************************************** * Residual Circulation ******************************************************** @@ -192,11 +190,13 @@ endif 'set t 1' 'd res'season'obs' -'set dfile ' zdifile +'set dfile 'zdifile 'set x 1' 'set t 1' 'sety' 'setz' +'getinfo undef' + undef = result if( CINTDIFF != NULL ) * -------------------- @@ -267,8 +267,19 @@ if( CINTDIFF != NULL ) else dn = 0 - 'shades 1' - 'd maskout( resdifz,abs(resdifz)-0.3 )' + 'define test = maskout( resdifz,abs(resdifz)-0.3 )' + 'minmax test' + testmax = subwrd(result,1) + testmin = subwrd(result,2) + if( testmax = undef & testmin = undef ) + 'shades resdifz 0' + cint = result + 'shades 'cint + 'd maskout( resdifz,abs(resdifz)-'cint' )' + else + 'shades 0.3' + 'd maskout( resdifz,abs(resdifz)-0.3 )' + endif endif @@ -278,7 +289,7 @@ endif 'set z 'z1' 'z2 'set t 1' -'cbarn -scale 0.9 -snum 0.5' +'cbarn -scale 0.9 -snum 0.5 -ymid 0.40' 'set gxout contour' 'set clopts 1 3 0.06' 'set ccolor 1' @@ -292,11 +303,14 @@ endif 'set parea off' 'set vpage 0 8.5 0.0 11' +'set string 1 l 4' +'set strsiz 0.07' +'draw string 0.05 0.10 ( EXPID: 'expid' )' + 'set string 1 c 6' 'set strsiz .11' -'draw string 4.25 10.85 EXPID: 'expid -'draw string 4.25 10.6 'expdsc' 'season' ('nmod')' +'draw string 4.25 10.75 'expdsc' 'season' ('nmod')' 'draw string 4.25 7.23 'obsdsc' 'season' ('nobs') ('climate')' if( dn != 0 ) diff --git a/plots/res/setup_epflx.gs b/plots/res/setup_epflx.gs index 89c81137..3a4e1403 100644 --- a/plots/res/setup_epflx.gs +++ b/plots/res/setup_epflx.gs @@ -3,11 +3,12 @@ function setup_epflx (args) source = subwrd(args,1) expid = subwrd(args,2) output = subwrd(args,3) +debug = subwrd(args,4) * Define Seasons to Process * ------------------------- seasons = '' - k = 4 + k = 5 while( k > 0 ) season = subwrd(args,k) if( season = '' ) @@ -24,7 +25,6 @@ say 'Running: setup_epflx 'source' 'expid' 'output' 'seasons say '--------------------------------------------------------' say ' ' - 'getinfo numfiles' numfiles = result @@ -43,42 +43,117 @@ if( numfiles = 'NULL' ) while( num <= numrc ) loc = num + 1 rcfile = subwrd( rcinfo,loc ) + say 'num = 'num' rcfile = 'rcfile '!grep filename 'rcfile' | cut -d'"\' -f2 > CTLFILE.txt" 'run getenv CTLFILE ' CTLFILE.num = result - '!echo `basename 'rcfile' | cut -d. -f2- > CMPID.txt`' - '!cat CMPID.txt | awk "{print length}" > LENID.txt' - 'run getenv LENID' - lenid = result - lenid = lenid - 3 - '!cat CMPID.txt | cut -b1-'lenid' > CMPID2.txt' - 'run getenv CMPID2' + '!basename 'rcfile' > BASENAME.txt' + 'run getenv BASENAME' + basename = result + length = strlen(basename) - 3 + say 'original length = 'length + '!echo 'basename' | cut -b1-'length' > CMPID.txt' + 'run getenv CMPID ' CMPID.num = result say ' CMPID #'num' = 'CMPID.num say 'CTLFILE #'num' = 'CTLFILE.num - 'open 'CTLFILE.num + + '!remove TEM_NAME.txt.' + '!basename 'CTLFILE.num' | cut -d. -f2 > TEM_NAME.txt' + say 'run getenv "TEM_NAME"' + 'run getenv "TEM_NAME"' + TEM_NAME.num = result + say 'TEM_Collection = 'TEM_NAME.num + pause + + 'xdfopen 'CTLFILE.num if( CMPID.num = expid ) ; nexpid = num ; endif num = num + 1 endwhile -say ' ' - + 'getinfo numfiles' numfiles = result endif -'run getenv MASKFILE ' maskfile - maskfile = result +* Loop over Experiment Datasets +* ----------------------------- +'getpwd' + pwd = result + +'getnumrc 'pwd + rcinfo = result + numrc = subwrd( rcinfo,1 ) + num = 1 + cnt = 0 +while( num <= numrc ) + loc = num + 1 + rcfile = subwrd( rcinfo,loc ) + say 'num = 'num' rcfile = 'rcfile + '!grep filename 'rcfile' | cut -d'"\' -f2 > CTLFILE.txt" + 'run getenv CTLFILE ' + CTLFILE.num = result + say 'CTLFILE.'num' = 'CTLFILE.num + + '!basename 'CTLFILE.num' > BASENAME.txt' + 'run getenv BASENAME' + basename = result + + '!remove NODE.txt' + '!echo 'basename' | cut -d. -f1 >> NODE.txt' + 'run getenv "NODE"' + node = result + cmpid = node + + say 'cmpid = 'cmpid + TEM_Collection = TEM_NAME.num + m = 2 + '!remove NODE.txt' + '!echo 'basename' | cut -d. -f'm' >> NODE.txt' + 'run getenv "NODE"' + node = result + while( node != TEM_Collection ) + cmpid = cmpid'.'node + say 'cmpid = 'cmpid + m = m + 1 + '!remove NODE.txt' + '!echo 'basename' | cut -d. -f'm' >> NODE.txt' + 'run getenv "NODE"' + node = result + endwhile + + CMPID.num = cmpid + say ' CMPID #'num' = 'CMPID.num + say 'CTLFILE #'num' = 'CTLFILE.num + say ' ' +* 'xdfopen 'CTLFILE.num + if( CMPID.num = expid ) ; nexpid = num ; endif + num = num + 1 +endwhile + +'getinfo numfiles' + numfiles = result + +'!/bin/rm -f LOCKFILE' 'q files' -say result + check = subwrd(result,1) -'run getenv "GEOSUTIL"' - geosutil = result +* Set MASKFILE to file which is not time-continuous (i.e., has UNDEF periods) +* --------------------------------------------------------------------------- + maskfile = 'NULL' +* maskfile = 2 + 'run setenv MASKFILE ' maskfile + 'run setenv NUMFILES ' numfiles -* Find time subset which includes all files -* ----------------------------------------- +say '' +say 'Files:' +'q files' +say result + +say 'Initializing BEGDATE and ENDATE for each Experiment:' +say '----------------------------------------------------' n = 1 while( n<=numfiles ) 'set dfile 'n @@ -92,7 +167,7 @@ while( n<=numfiles ) 'getinfo date' enddate.n = result 'run setenv ENDDATE'n' 'enddate.n -say 'begdate.'n': 'begdate.n' enddate.'n': 'enddate.n + say 'begdate.'n': 'begdate.n' enddate.'n': 'enddate.n n = n + 1 endwhile @@ -123,6 +198,8 @@ while( n<=numfiles ) n = n + 1 endwhile +* Hardwire begdate and enddate here +* --------------------------------- *begdateA = 00z01jan2010 *enddateA = 00z01dec2014 @@ -144,6 +221,8 @@ say 'set time 'begdateA' 'enddateA 'set time 'begdateA' 'enddateA say ' ' +say 'setdates' + 'setdates' say 'getdates' 'getdates' say ' ' @@ -165,6 +244,7 @@ while( n<=numfiles ) 'seasonal epfy A'n 'seasonal epfz A'n 'seasonal epfdiv A'n + n = n + 1 endwhile @@ -176,6 +256,9 @@ endwhile 'set gxout shaded' 'c' +'run getenv "GEOSUTIL"' + geosutil = result + n = 1 while( n<=numfiles ) k = 1 @@ -185,9 +268,22 @@ while( n<=numfiles ) k = k+1 say 'Running: epflx.gs 'CMPID.n' 'season' A'n' 'output say '-------------------------------------------------' - 'run 'geosutil'/plots/res/epflx.gs 'CMPID.n' 'season' A'n' 'output - pause - 'c' + 'set x 1' + 'sety' + 'setz' + flag = "" + while ( flag = "" ) + 'run 'geosutil'/plots/res/epflx.gs 'CMPID.n' 'TEM_NAME.n' 'season' A'n' 'output + if( debug = "debug" ) + say "Hit ENTER to repeat plot" + say "Type 'next' for next plot, 'done' for next field" + pull flag + else + flag = "next" + endif + 'c' + endwhile + else k = -1 endif @@ -206,14 +302,35 @@ if( CMPID.n != expid ) k = k+1 say 'Running: epflx_diff.gs 'expid' 'CMPID.n' 'season' A'nexpid' A'n' 'output say '----------------------------------------------------' - 'run 'geosutil'/plots/res/epflx_diff.gs 'expid' 'CMPID.n' 'season' A'nexpid' A'n' 'output - pause - 'c' - say 'Running: epflx_diff.gs 'CMPID.n' 'expid' 'season' A'n' A'nexpid' 'output - say '----------------------------------------------------' - 'run 'geosutil'/plots/res/epflx_diff.gs 'CMPID.n' 'expid' 'season' A'n' A'nexpid' 'output - pause - 'c' + + flag = "" + while ( flag = "" ) + say 'run 'geosutil'/plots/res/epflx_diff.gs 'expid' 'CMPID.n' 'TEM_NAME.nexpid' 'TEM_NAME.n' 'season' A'nexpid' A'n' 'output + 'run 'geosutil'/plots/res/epflx_diff.gs 'expid' 'CMPID.n' 'TEM_NAME.nexpid' 'TEM_NAME.n' 'season' A'nexpid' A'n' 'output + if( debug = "debug" ) + say "Hit ENTER to repeat plot" + say "Type 'next' for next plot, 'done' for next field" + pull flag + else + flag = "next" + endif + 'c' + endwhile + + flag = "" + while ( flag = "" ) + say 'run 'geosutil'/plots/res/epflx_diff.gs 'CMPID.n' 'expid' 'TEM_NAME.n' 'TEM_NAME.nexpid' 'season' A'n' A'nexpid' 'output + 'run 'geosutil'/plots/res/epflx_diff.gs 'CMPID.n' 'expid' 'TEM_NAME.n' 'TEM_NAME.nexpid' 'season' A'n' A'nexpid' 'output + if( debug = "debug" ) + say "Hit ENTER to repeat plot" + say "Type 'next' for next plot, 'done' for next field" + pull flag + else + flag = "next" + endif + 'c' + endwhile + else k = -1 endif @@ -224,10 +341,10 @@ endwhile * Move DATA for Archive * --------------------- -'!/bin/mv -f residual*ctl ../' -'!/bin/mv -f residual*data ../' -'!/bin/mv -f VERIFICATION*rc ../' -'!/bin/cp -f ../residual.'expid'.* ../../' +*'!/bin/mv -f residual*ctl ../' +*'!/bin/mv -f residual*data ../' +*'!/bin/mv -f VERIFICATION*rc ../' +*'!/bin/cp -f ../residual.'expid'.* ../../' return diff --git a/plots/res/setup_wstar.gs b/plots/res/setup_wstar.gs index 9905c1a3..92c5951a 100644 --- a/plots/res/setup_wstar.gs +++ b/plots/res/setup_wstar.gs @@ -39,7 +39,6 @@ if( numfiles = 'NULL' ) * Loop over Experiment Datasets * ----------------------------- -say '' 'getpwd' pwd = result @@ -67,7 +66,16 @@ while( num <= numrc ) say ' CMPID #'num' = 'CMPID.num say 'CTLFILE #'num' = 'CTLFILE.num - 'open 'CTLFILE.num + + '!remove TEM_NAME.txt.' + '!basename 'CTLFILE.num' | cut -d. -f2 > TEM_NAME.txt' + say 'run getenv "TEM_NAME"' + 'run getenv "TEM_NAME"' + TEM_NAME.num = result + say 'TEM_Collection = 'TEM_NAME.num + pause + + 'xdfopen 'CTLFILE.num if( CMPID.num = expid ) ; nexpid = num ; endif num = num + 1 endwhile @@ -89,9 +97,6 @@ endif 'run setenv NUMFILES ' numfiles say '' -say 'Files:' -'q files' -say result say 'MASKFILE: 'maskfile say '' @@ -223,6 +228,7 @@ while( n<=numfiles ) if( maskfile != 'NULL' ) 'define wstar'n' = maskout( wstar.'n',abs(wstar.'maskfile') )' + 'define res'n' = maskout( res.'n' ,abs( res.'maskfile') )' 'define wstar = maskout( wstar.'n',abs(wstar.'maskfile') )' 'define res = maskout( res.'n' ,abs( res.'maskfile') )' * 'define epfy = maskout( epfy.'n' ,abs( epfy.'maskfile') )' @@ -232,10 +238,12 @@ while( n<=numfiles ) * 'define vstar = maskout( vstar.'n',abs(vstar.'maskfile') )' else 'define wstar'n' = wstar.'n + 'define res'n' = res.'n endif 'seasonal wstar A'n 'seasonal res A'n + * 'seasonal epfy A'n * 'seasonal epfz A'n * 'seasonal epfdiv A'n @@ -317,11 +325,14 @@ while( k > 0 ) n = 1 while( n<=numfiles ) + CMPID = CMPID.n + TEM_Collection = TEM_NAME.n + 'set dfile 'n 'getinfo desc' desc = result node = '' - m = 2 + m = 1 '!remove NODE.txt' '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' @@ -333,7 +344,7 @@ while( k > 0 ) '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' node = result - while( node != 'ctl' & node != 'data' ) + while( node != TEM_Collection ) EXP.n = EXP.n'.'node say 'EXP'n' = 'EXP.n m = m + 1 diff --git a/plots/res/turn_around_lats.gs b/plots/res/turn_around_lats.gs index 36b045d3..d387f574 100644 --- a/plots/res/turn_around_lats.gs +++ b/plots/res/turn_around_lats.gs @@ -34,6 +34,51 @@ levmin = 20 'run getenv TIME ' time = result + +* Loop over Experiment Datasets to get Experiment IDs +* --------------------------------------------------- +say '' +'getpwd' + pwd = result + +'getnumrc 'pwd + rcinfo = result + numrc = subwrd( rcinfo,1 ) + num = 1 + cnt = 0 +while( num <= numrc ) + loc = num + 1 + rcfile = subwrd( rcinfo,loc ) + say 'num = 'num' rcfile = 'rcfile + + '!grep filename 'rcfile' | cut -d'"\' -f2 > CTLFILE.txt" + 'run getenv CTLFILE ' + CTLFILE.num = result + + '!basename 'rcfile' > BASENAME.txt' + 'run getenv BASENAME' + basename = result + length = strlen(basename) - 3 + say 'original length = 'length + '!echo 'basename' | cut -b1-'length' | cut -d. -f2 > CMPID.txt' + 'run getenv CMPID ' + CMPID.num = result + + say ' CMPID #'num' = 'CMPID.num + say ' CTLFILE #'num' = 'CTLFILE.num + + '!remove TEM_NAME.txt.' + '!basename 'CTLFILE.num' | cut -d. -f2 > TEM_NAME.txt' + say 'run getenv "TEM_NAME"' + 'run getenv "TEM_NAME"' + TEM_NAME.num = result + say 'TEM_Collection = 'TEM_NAME.num + pause + + num = num + 1 +endwhile + + 'getinfo numfiles' numfiles = result @@ -63,6 +108,7 @@ levmin = 20 'set datawarn off' nbeg = 1 + nmax = numfiles 'run getenv BEGDATE' @@ -82,6 +128,8 @@ filecount = 0 * Loop over Files * --------------- +'q files' +say 'FILES: 'result say 'Looping over files 'nbeg' to 'nmax' to find TA Lats' n =nbeg while( n<=nmax ) @@ -245,7 +293,7 @@ n = n + 1 endwhile say 'Average wstr'season'latmn = ' wstrlatmn ' wstr'season'latmx = ' wstrlatmx say ' ' -*pause + pause * ---------------------------------------------------- @@ -284,6 +332,7 @@ ybot_wstar = subwrd(result,2) 'getinfo dlat' dlat1 = result +found = TRUE count = 0 n =nbeg while( n<=nmax ) @@ -296,6 +345,8 @@ if( dlat = dlat1 & foundS.n = TRUE & foundN.n = TRUE ) 'set dfile 1' 'define wstrave = wstrave + 1000*wstr'season''time''n count = count + 1 +else + found = FALSE endif n = n + 1 endwhile @@ -326,21 +377,29 @@ while( n<=nmax ) 'set cmark 0' 'set cthick 1' 'd 1000*wstr'season''time''n +say 'Turn-Around Lats for File 'n' Color: 'color.n +pause n = n + 1 endwhile 'set dfile 1' +if( found = TRUE ) 'set cstyle 1' 'set cmark 0' 'set cthick 8' 'set ccolor 1' 'd wstrave' +say 'Turn-Around Lats for Average Color: '1 +pause +endif 'set cmark 0' 'set ccolor 1' 'set cthick 1' 'set z 1' 'd lev-lev' +say 'Zero Line Color: '1 +pause * Compute Turn-Around Lats based on wstrave * ----------------------------------------- @@ -435,12 +494,14 @@ while( n<=nmax ) n = n + 1 endwhile +if( found = TRUE ) 'set dfile 1' 'set cstyle 1' 'set cmark 0' 'set cthick 8' 'set ccolor 1' 'd psiave' +endif 'set cmark 0' 'set ccolor 1' @@ -537,13 +598,19 @@ ymax = subwrd(result,6) '!remove DESC.txt' '!echo 'numfiles' > turn_around_'season'.stack' +'q file' +say 'FILES: 'result +pause n = 1 while( n<=numfiles ) + CMPID = CMPID.n + TEM_Collection = TEM_NAME.n + 'set dfile 'n 'getinfo desc' desc = result node = '' - m = 2 + m = 1 '!remove NODE.txt' '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' @@ -555,7 +622,11 @@ ymax = subwrd(result,6) '!basename 'desc' | cut -d. -f'm' >> NODE.txt' 'run getenv "NODE"' node = result - while( node != 'ctl' & node != 'data' ) + + say 'NODE = 'node + say 'TEM_Collection = 'TEM_Collection + + while( node != TEM_Collection ) EXP.n = EXP.n'.'node say 'EXP'n' = 'EXP.n m = m + 1 diff --git a/plots/res/zonal.gs b/plots/res/zonal.gs index 6f7f1e60..253ed6cc 100644 --- a/plots/res/zonal.gs +++ b/plots/res/zonal.gs @@ -22,15 +22,16 @@ function zonal (args) say ' EXPID: 'expid say 'EXPDSC: 'expdsc +say ' ' 'run getenv "GEOSUTIL"' - geosutil = result + geosutil = result 'run getenv "VERIFICATION"' - verification = result + verification = result 'run getenv "ARCH"' - arch = result + arch = result 'run getpwd' pwd = result @@ -38,105 +39,110 @@ say 'EXPDSC: 'expdsc '!/bin/cp 'geosutil'/plots/res/VERIFICATION*rc .' * ------------------------------------------------------------------------------------------------- -* Note: Assume that we ALWAYS want to create an updated residual.EXPID.ctl and residual.EXPID.data +* Check for TEM_Collection Diagnostic Files under source home directory * ------------------------------------------------------------------------------------------------- -* Check for residual.EXPID.data under source home directory -* --------------------------------------------------------- 'run getenv "SOURCE"' source = result -say 'Checking for 'source'/residual.'expid'.ctl' + +'run getenv "HISTORYRC"' + HISTORYRC = result + +'!remove TEM_NAME' +'!/bin/grep TEM_Collection 'HISTORYRC' > TEM_NAME' + dummy = read(TEM_NAME) + ioflag = sublin( dummy,1 ) + record = sublin( dummy,2 ) + say expid' ioflag: 'ioflag + say expid' record: 'record + if( ioflag = 0 ) + TEM_Collection = subwrd( record,2 ) + else + TEM_Collection = 'TEM_Diag' + endif + dummy = close(TEM_NAME) +say ' ' +say 'EXPID Transformed Eulerian Mean (TEM) Diagnostic Collection: 'TEM_Collection +say '----------------------------------------------------------- ' +say ' ' +say 'Checking for 'source'/'TEM_Collection'/'expid'.'TEM_Collection'.monthly.ddf' '!remove CHECKFILE.txt' -'!chckfile 'source'/residual.'expid'.ctl' - 'run getenv CHECKFILE' - expid.ctl = result -'!remove CHECKFILE.txt' -'!chckfile 'source'/residual.'expid'.data' +'!chckfile 'source'/'TEM_Collection'/'expid'.'TEM_Collection'.monthly.ddf' 'run getenv CHECKFILE' - expid.data = result - expid.ctl = NULL - pause EXPID.CTL = expid.ctl + expid.ddf = result + say 'expid.ddf: 'expid.ddf + pause -* Check for residual.EXPID.data under plot Output directory -* --------------------------------------------------------- -if( expid.ctl = 'NULL' ) - say 'Checking for 'pwd'/../residual.'expid'.ctl' +if( expid.ddf = 'NULL' ) +* Check for old format residual.EXPID.data under plot Output directory +* -------------------------------------------------------------------- + say 'Checking for 'pwd'/../residual.'expid'.ctl' '!remove CHECKFILE.txt' '!chckfile 'pwd'/../residual.'expid'.ctl' 'run getenv CHECKFILE' expid.ctl = result + '!remove CHECKFILE.txt' '!chckfile 'pwd'/../residual.'expid'.data' 'run getenv CHECKFILE' expid.data = result - pause EXPID.DATA = expid.data -endif - -* Check EP_UPDATE to over-ride defaults -* ------------------------------------- - 'run getenv EP_UPDATE' - EP_UPDATE = result -* Setting EP_UPDATE = TRUE to force re-calculation -* ------------------------------------------------ -* EP_UPDATE = TRUE + if( ( expid.ctl = 'NULL' & expid.data != 'NULL' ) | ( expid.ctl != 'NULL' & expid.data = 'NULL' ) ) + say 'Inconsistent 'expid' ctl and data files:' + say 'expid.ctl : 'expid.ctl + say 'expid.data: 'expid.data + return + endif -if( expid.ctl = 'NULL' | EP_UPDATE != 'NULL' ) + if( expid.ctl != 'NULL' & expid.data != 'NULL' ) + '!/bin/cp 'pwd'/../residual.'expid'.ctl .' + '!/bin/cp 'pwd'/../residual.'expid'.data .' -* Updated method to overwrite existing residual.EXPID.data -* -------------------------------------------------------- - say '!/bin/mv -f 'source'/residual.'expid'.ctl 'source'/residual.'expid'.ctl.old ' - say '!/bin/mv -f 'source'/residual.'expid'.data 'source'/residual.'expid'.data.old ' +* Create TEM_Collection Diagnostic data based on old format +* --------------------------------------------------------- + '!'geosutil'/plots/grads_util/Create_TEM_ddf.csh 'expid' 'TEM_Collection - '!remove CHECKFILE.txt' - '!chckfile 'source'/residual.'expid'.ctl' - 'run getenv CHECKFILE' - CHECKFILE = result - if( CHECKFILE != 'NULL' ) - '!/bin/mv -f 'source'/residual.'expid'.ctl 'source'/residual.'expid'.ctl.old ' - endif + say 'Create TEM_Collection Directory under SOURCE' +* ------------------------------------------------- + '!/bin/mkdir -p 'source'/'TEM_Collection + '!/bin/cp 'expid'.'TEM_Collection'.monthly.* 'source'/'TEM_Collection - '!remove CHECKFILE.txt' - '!chckfile 'source'/residual.'expid'.data' - 'run getenv CHECKFILE' - CHECKFILE = result - if( CHECKFILE != 'NULL' ) - '!/bin/mv -f 'source'/residual.'expid'.data 'source'/residual.'expid'.data.old ' - endif + say 'Move TEM_Collection Diagnostic data to OUTPUT directory' +* ------------------------------------------------------------ + '!/bin/mv 'expid'.'TEM_Collection'.monthly.* ../' - expid.ctl = 'NULL' + '!remove CHECKFILE.txt' + '!chckfile ../'expid'.'TEM_Collection'.monthly.ddf' + 'run getenv CHECKFILE' + expid.ddf = result + endif + pause +endif + +* ------------------------------------------------------------------------------------------------- +* Using existing EXPID.TEM_Collection data to create corresponding VERIFICATION rc file +* ------------------------------------------------------------------------------------------------- -else -* Using existing residual.EXPID.data -* ---------------------------------- - say 'Using existing residual.EXPID.data' -* '!/bin/cp 'source'/residual.'expid'.ctl .' -* '!/bin/cp 'source'/residual.'expid'.data .' - '!/bin/cp 'expid.ctl' .' - '!/bin/cp 'expid.data' .' +if( expid.ddf != 'NULL' ) + say 'Creating VERIFICATION.expid.rc associated with existing 'expid'.'TEM_Collection'.monthly.ddf' '!remove sedfile' '!touch sedfile' - '!echo "s?@EXPID?"'expid'?g >> sedfile' - '!echo "s?@EXPDSC?"'expdsc'?g >> sedfile' - '!echo "s?@pwd?"'pwd'?g >> sedfile' + '!echo "s?@EXPID?"'expid'?g >> sedfile' + '!echo "s?@EXPDSC?"'expdsc'?g >> sedfile' + '!echo "s?@filename?"'expid.ddf'?g >> sedfile' '!/bin/cp 'geosutil'/plots/res/VERIFICATION.rc.tmpl .' '!sed -f sedfile VERIFICATION.rc.tmpl > VERIFICATION.'expid'.rc' '!remove VERIFICATION.rc.tmpl' - endif + pause +endif -say ' ' -say 'expid: 'expid -say 'expid.ctl: 'expid.ctl -say 'expid.data: 'expid.data -say ' ' -*pause +* ------------------------------------------------------------------------------------------------- +* Compute EXPID TEM_Diagnostic calculations if needed + if( expid.ddf = 'NULL' ) +* ------------------------------------------------------------------------------------------------- -* Perform residual calculations -* ----------------------------- -if( expid.ctl = 'NULL' ) - * Initialize Environment using V-Wind File * ---------------------------------------- 'set dfile 'vfile @@ -357,30 +363,52 @@ t=t+1 endwhile 'disable fwrite' -* Run Fortran Code to Produce StreamFunction and Residual Circulation -* ------------------------------------------------------------------- -say ' 'geosutil'/bin/zonal_'arch'.x -tag 'expid' -desc 'descm -'! 'geosutil'/bin/zonal_'arch'.x -tag 'expid' -desc 'descm - -'!remove sedfile' -'!touch sedfile' - -'!echo "s?@EXPID?"'expid'?g >> sedfile' -'!echo "s?@EXPDSC?"'expdsc'?g >> sedfile' -'!echo "s?@pwd?"'pwd'?g >> sedfile' -'!/bin/cp 'geosutil'/plots/res/VERIFICATION.rc.tmpl .' -'!sed -f sedfile VERIFICATION.rc.tmpl > VERIFICATION.'expid'.rc' -'!remove VERIFICATION.rc.tmpl' - -* Copy DATA to plot Output and Source Directory -* --------------------------------------------- -'!/bin/cp -f residual.'expid'.ctl ../' -'!/bin/cp -f residual.'expid'.data ../' -'!/bin/cp -f VERIFICATION.'expid'.rc ../' -'!/bin/cp -f ../residual.'expid'.* ../../' +* Run Fortran Code to Produce Old Format residual.EXPID.ctl and residual.EXPID.data +* --------------------------------------------------------------------------------- +say ' 'geosutil'/plots/zonal_'arch'.x -tag 'expid' -desc 'descm +'! 'geosutil'/plots/zonal_'arch'.x -tag 'expid' -desc 'descm -endif -* ----------------------------------------------------- + +* Create TEM_Collection Diagnostic data based on old format +* --------------------------------------------------------- + '!'geosutil'/plots/grads_util/Create_TEM_ddf.csh 'expid' 'TEM_Collection + + +say 'Create TEM_Collection Directory under SOURCE' +* ----------------------------------------------- + '!/bin/mkdir -p 'source'/'TEM_Collection + '!/bin/cp 'expid'.'TEM_Collection'.monthly.* 'source'/'TEM_Collection + +say 'Move TEM_Collection Diagnostic data to OUTPUT directory' +* ---------------------------------------------------------- + '!/bin/mv 'expid'.'TEM_Collection'.monthly.* ../' + + '!remove CHECKFILE.txt' + '!chckfile ../'expid'.'TEM_Collection'.monthly.ddf' + 'run getenv CHECKFILE' + expid.ddf = result + +say 'Creating VERIFICATION.expid.rc associated with 'expid'.'TEM_Collection'.monthly.ddf' + '!remove sedfile' + '!touch sedfile' + '!echo "s?@EXPID?"'expid'?g >> sedfile' + '!echo "s?@EXPDSC?"'expdsc'?g >> sedfile' + '!echo "s?@filename?"'expid.ddf'?g >> sedfile' + '!/bin/cp 'geosutil'/plots/res/VERIFICATION.rc.tmpl .' + '!sed -f sedfile VERIFICATION.rc.tmpl > VERIFICATION.'expid'.rc' + '!remove VERIFICATION.rc.tmpl' + pause + +* ------------------------------------------------------------------------------------------------- +* End Computation of EXPID TEM_Diagnostics + endif +* ------------------------------------------------------------------------------------------------- + + + +* ------------------------------------------------------------------------------------------------- +* Loop over Possible Experiment Datasets for Comparison +* ------------------------------------------------------------------------------------------------- 'run setdates' @@ -400,6 +428,13 @@ say ' ' say 'Comparing with: 'exp say 'Comparison type: 'type + '!remove NODE.txt' + '!basename 'exp' >> NODE.txt' + 'run getenv "NODE"' + cmpid = result +say ' CMPID: 'cmpid +pause + '!remove CHECKFILE.txt' '!chckfile 'exp'/.HOMDIR' 'run getenv CHECKFILE' @@ -410,6 +445,25 @@ say 'Comparison type: 'type '!/bin/cp 'exp'/HISTORY.rc .' endif +'!remove TEM_NAME' +'!/bin/grep TEM_Collection HISTORY.rc > TEM_NAME' + dummy = read(TEM_NAME) + ioflag = sublin( dummy,1 ) + record = sublin( dummy,2 ) + say cmpid' ioflag: 'ioflag + say cmpid' record: 'record + if( ioflag = 0 ) + TEM_Collection = subwrd( record,2 ) + else + TEM_Collection = 'TEM_Diag' + endif + dummy = close(TEM_NAME) +say ' ' +say 'CMPID Transformed Eulerian Mean (TEM) Diagnostic Collection: 'TEM_Collection +say '----------------------------------------------------------- ' +say ' ' +pause + '!cat HISTORY.rc | sed -e "s/,/ , /g" | sed -e "s/*/@/g" > HISTORY.T' 'run getvar V DYN 'exp @@ -434,59 +488,89 @@ say 'Comparison Desc: 'obsdsc pause -* Check for residual.OBSID.data in comparison experiment directory -* Note: if residual.OBSID.data exists, then create an associated VERIFICATION.obsid.rc -* ------------------------------------------------------------------------------------- - say 'Checking for 'exp'/residual.'obsid'.ctl' -'!remove CHECKFILE.txt' -'!chckfile 'exp'/residual.'obsid'.ctl' - 'run getenv CHECKFILE' - obsid.ctl = result +* Check for Transport Diagnostic Files under comparison experiment directory +* -------------------------------------------------------------------------- + say 'Checking for 'exp'/'TEM_Collection'/'obsid'.'TEM_Collection'.monthly.ddf' + '!remove CHECKFILE.txt' -'!chckfile 'exp'/residual.'obsid'.data' +'!chckfile 'exp'/'TEM_Collection'/'obsid'.'TEM_Collection'.monthly.ddf' 'run getenv CHECKFILE' - obsid.data = result - obsid.ctl = NULL + obsid.ddf = result + say 'obsid.ddf: 'obsid.ddf +pause -* Check for residual.OBSID.data in plot Output directory -* ------------------------------------------------------ -if( obsid.ctl = 'NULL' ) +* Check for old format residual.OBSID.data in plot Output directory +* ----------------------------------------------------------------- +if( obsid.ddf = 'NULL' ) say 'Checking for 'pwd'/../residual.'obsid'.ctl' + '!remove CHECKFILE.txt' '!chckfile 'pwd'/../residual.'obsid'.ctl' 'run getenv CHECKFILE' obsid.ctl = result + '!remove CHECKFILE.txt' '!chckfile 'pwd'/../residual.'obsid'.data' 'run getenv CHECKFILE' obsid.data = result + + if( ( obsid.ctl = 'NULL' & obsid.data != 'NULL' ) | ( obsid.ctl != 'NULL' & obsid.data = 'NULL' ) ) + say 'Inconsistent 'obsid' ctl and data files:' + say 'expid.ctl : 'expid.ctl + say 'expid.data: 'expid.data + return + endif + + if( obsid.ctl != 'NULL' & obsid.data != 'NULL' ) + '!/bin/cp 'pwd'/../residual.'obsid'.ctl .' + '!/bin/cp 'pwd'/../residual.'obsid'.data .' + + '!'geosutil'/plots/grads_util/Create_TEM_ddf.csh 'obsid' 'TEM_Collection + + '!remove CHECKFILE.txt' + '!chckfile 'obsid'.'TEM_Collection'.monthly.ddf' + 'run getenv CHECKFILE' + obsid.ddf = result + say 'obsid.ddf: 'obsid.ddf + pause + endif endif -* Setting obsid.ctl = NULL to force re-calculation -* ------------------------------------------------ -* obsid.ctl = 'NULL' +* Check for new format 'obsid'.'TEM_Collection'.monthly.ddf in plot Output directory +* ---------------------------------------------------------------------------------- +if( obsid.ddf = 'NULL' ) + say 'Checking for 'pwd'/../'obsid'.'TEM_Collection'.monthly.ddf' - if( obsid.ctl != 'NULL' ) + '!remove CHECKFILE.txt' + '!chckfile 'pwd'/../'obsid'.'TEM_Collection'.monthly.ddf' + 'run getenv CHECKFILE' + obsid.ddf = result + say 'obsid.ddf: 'obsid.ddf + pause +endif - say 'Creating VERIFICATION.'obsid'.rc for 'pwd'/../residual.'obsid'.data' -* '!/bin/cp 'pwd'/../residual.'obsid'.ctl .' -* '!/bin/cp 'pwd'/../residual.'obsid'.data .' - '!/bin/cp 'obsid.ctl' .' - '!/bin/cp 'obsid.data' .' +* ------------------------------------------------------------------------------------- +* Create an associated VERIFICATION.obsid.rc +* ------------------------------------------------------------------------------------- + + if( obsid.ddf != 'NULL' ) + say 'Creating VERIFICATION.obsid.rc associated with 'obsid'.'TEM_Collection'.monthly.ddf' '!remove sedfile' '!touch sedfile' - '!echo "s?@EXPID?"'obsid'?g >> sedfile' - '!echo "s?@EXPDSC?"'obsdsc'?g >> sedfile' - '!echo "s?@pwd?"'pwd'?g >> sedfile' + '!echo "s?@EXPID?"'obsid'?g >> sedfile' + '!echo "s?@EXPDSC?"'obsdsc'?g >> sedfile' + '!echo "s?@filename?"'obsid.ddf'?g >> sedfile' '!/bin/cp 'geosutil'/plots/res/VERIFICATION.rc.tmpl .' '!sed -f sedfile VERIFICATION.rc.tmpl > VERIFICATION.'obsid'.rc' '!remove VERIFICATION.rc.tmpl' - - endif + endif + +* ------------------------------------------------------------------------------------- +* ------------------------------------------------------------------------------------- * Checking for VERIFICATION.obsid.rc -* Note: if VERIFICATION.obsid.rc exists, then its contents should point to residual.'obsid'.ctl -* --------------------------------------------------------------------------------------------- +* Note: if VERIFICATION.obsid.rc exists, then its contents should point to 'obsid'.'TEM_Collection'.monthly.ddf +* ---------------------------------------------------------------------------------------------------------- '!remove CHECKFILE.txt' '!chckfile VERIFICATION.'obsid'.rc' 'run getenv CHECKFILE' @@ -497,21 +581,19 @@ endif endif -* Setting obsid.rc = NULL to force re-calculation -* ----------------------------------------------- -* obsid.rc = 'NULL' - say ' ' say 'obsid: 'obsid -say 'obsid.ctl: 'obsid.ctl -say 'obsid.data: 'obsid.data say 'obsid.rc: 'obsid.rc say ' ' -*pause +pause + +* ------------------------------------------------------------------------------------------------- +* Compute OBSID TEM_Diagnostic calculations if needed + if( obsid != 'ERA5' & obsid.ddf = 'NULL' ) +* ------------------------------------------------------------------------------------------------- * Perform residual calculations * ----------------------------- - if( obsid != 'ERA5' & ( obsid.ctl = 'NULL' | obsid.rc = 'NULL' ) ) 'set dfile 'vfile 'set y 1' @@ -715,27 +797,45 @@ t=t+1 endwhile 'disable fwrite' -* Run Fortran Code to Produce StreamFunction and Residual Circulation -* ------------------------------------------------------------------- -say ' 'geosutil'/bin/zonal_'arch'.x -tag 'obsid' -desc 'desco -'! 'geosutil'/bin/zonal_'arch'.x -tag 'obsid' -desc 'desco -'!remove sedfile' -'!touch sedfile' +'!pwd' +say 'Current Direcory: 'result +pause -'!echo "s?@EXPID?"'obsid'?g >> sedfile' -'!echo "s?@EXPDSC?"'obsdsc'?g >> sedfile' -'!echo "s?@pwd?"'pwd'?g >> sedfile' -'!/bin/cp 'geosutil'/plots/res/VERIFICATION.rc.tmpl .' -'!sed -f sedfile VERIFICATION.rc.tmpl > VERIFICATION.'obsid'.rc' -'!remove VERIFICATION.rc.tmpl' -'!cat sedfile' +* Run Fortran Code to Produce Old Format residual.OBSID.ctl and residual.OBSID.data +* --------------------------------------------------------------------------------- +say ' 'geosutil'/plots/zonal_'arch'.x -tag 'obsid' -desc 'desco +'! 'geosutil'/plots/zonal_'arch'.x -tag 'obsid' -desc 'desco + +pause + +* Create TEM_Collection Diagnostic data based on old format +* --------------------------------------------------------- + '!'geosutil'/plots/grads_util/Create_TEM_ddf.csh 'obsid' 'TEM_Collection + + '!remove CHECKFILE.txt' + '!chckfile 'obsid'.'TEM_Collection'.monthly.ddf' + 'run getenv CHECKFILE' + obsid.ddf = result +pause + +say 'Creating VERIFICATION.obsid.rc associated with 'obsid'.'TEM_Collection'.monthly.ddf' +* -------------------------------------------------------------------------------- + '!remove sedfile' + '!touch sedfile' + '!echo "s?@EXPID?"'obsid'?g >> sedfile' + '!echo "s?@EXPDSC?"'obsdsc'?g >> sedfile' + '!echo "s?@filename?"'obsid.ddf'?g >> sedfile' + '!/bin/cp 'geosutil'/plots/res/VERIFICATION.rc.tmpl .' + '!sed -f sedfile VERIFICATION.rc.tmpl > VERIFICATION.'obsid'.rc' + '!remove VERIFICATION.rc.tmpl' +pause + +say 'Copying 'obsid' TEM_Collection Diagnostic data to PLOT directory' +* ------------------------------------------------------------------- + '!/bin/cp 'obsid'.'TEM_Collection'.monthly.* ../' +pause -* Copy DATA to Output Directory -* ----------------------------- -'!/bin/cp -f residual.'obsid'.ctl ../' -'!/bin/cp -f residual.'obsid'.data ../' -'!/bin/cp -f VERIFICATION.'obsid'.rc ../' * End VERIFICATION.obsid.rc CHECKFILE Test * ---------------------------------------- diff --git a/plots/zcmp/zcloseplt b/plots/zcmp/zcloseplt index 72bd7c79..77c77593 100644 --- a/plots/zcmp/zcloseplt +++ b/plots/zcmp/zcloseplt @@ -143,8 +143,7 @@ if( scale = LOG ) else zid = 'z'ptop endif - -dpct = 'NULL' + dpct = 0.1 if( EXPORT = "U" ) oname = "/horiz_"cmpnam"_uwnd_"zid"_closeness_"obsnam @@ -271,6 +270,11 @@ if( EXPORT = "O3" ) 'getresource 'PLOTRC' O3_CHEMISTRY_FIXED_PLOT_FACTOR' ; fixpltfact = result 'getresource 'PLOTRC' O3_CHEMISTRY_FIXED_PLOT_CINT' ; fixpltcint = result endif + +if( dpct = 'NULL' ) + dpct = 0.1 +endif + say ' TITLE: 'title say ' CLEVS: 'clevs say ' DLEVS: 'dlevs @@ -499,12 +503,8 @@ if( scale = LOG ) ; 'setlevs' ; endif dqrel = result/100 say 'int dqrel = 'dqrel - if( dpct = 'NULL' ) - dpct = 0.1 - else - 'd 'dpct - dpct = subwrd(result,4) - endif + 'd 'dpct + dpct = subwrd(result,4) say 'DQMAX: 'qmax' QMOD_MAX: 'amdifmax say 'Relative Percent Difference: 'dqrel' (100*DQMAX/QMOD_MAX)' diff --git a/plots/zcmp/zplt b/plots/zcmp/zplt index 0fd2c649..db773bba 100644 --- a/plots/zcmp/zplt +++ b/plots/zcmp/zplt @@ -306,7 +306,7 @@ if( Lbit = L ) ; 'setlevs' ; endif 'set clevs 'clevs 'set ccols 'ccols 'd modz' -'draw ylab Pressure (mb)' +'draw ylab Pressure (hPa)' 'set gxout contour' 'set ccolor 1' 'set clevs 'clevs @@ -334,7 +334,7 @@ if( Lbit = L ) ; 'setlevs' ; endif 'set clevs 'clevs 'set ccols 'ccols 'd obsz' -'draw ylab Pressure (mb)' +'draw ylab Pressure (hPa)' 'set gxout contour' 'set clevs 'clevs 'set ccolor 1' @@ -420,7 +420,7 @@ endif 'set t 1' if( CINTDIFF = 'NULL' ) - 'draw ylab Pressure (mb)' + 'draw ylab Pressure (hPa)' endif 'set vpage off' @@ -428,7 +428,7 @@ endif if( CINTDIFF != 'NULL' ) 'set string 1 l 6 90' 'set strsiz .16' - 'draw string 0.75 1.47 Pressure (mb)' + 'draw string 0.75 1.47 Pressure (hPa)' endif 'set string 1 l 4 0' diff --git a/post/CMakeLists.txt b/post/CMakeLists.txt index dca4991a..370875e2 100644 --- a/post/CMakeLists.txt +++ b/post/CMakeLists.txt @@ -19,7 +19,7 @@ target_compile_definitions(post_nompi PRIVATE mpi) foreach (basename makeiau hdf2rs gg2fv eta2prs eta2rst rs_hinterp - time_ave convert_eta rs2hdf + convert_eta rs2hdf ec_prs2fv flat2hdf fvrst gg2eta ec_prs2eta merra2scm rsg3_vinterp era5_prs2eta) ecbuild_add_executable (TARGET ${basename}.x SOURCES ${basename}.F DEFINITIONS mpi) @@ -28,6 +28,8 @@ ecbuild_add_executable (TARGET binarytile.x SOURCES binarytile.F90) ecbuild_add_executable (TARGET rs_numtiles.x SOURCES rs_numtiles.F90) ecbuild_add_executable (TARGET rs_vinterp.x SOURCES rs_vinterp.F90) ecbuild_add_executable (TARGET rs_vinterp_scm.x SOURCES rs_vinterp_scm.F90) +ecbuild_add_executable (TARGET time_ave.x SOURCES time_ave.F TEM.F90 DEFINITIONS mpi) +ecbuild_add_executable (TARGET Create_TEM_Diag.x SOURCES Create_TEM_Diag.F TEM.F90) ecbuild_add_executable (TARGET stats.x SOURCES stats.F90) if (DISABLE_FIELD_WIDTH_WARNING) @@ -77,6 +79,7 @@ set_property(SOURCE merra2scm.F APPEND_STRING PROPERTY COMPILE_FLAGS "${BYTERECL set_property(SOURCE rs2hdf.F APPEND_STRING PROPERTY COMPILE_FLAGS "${BYTERECLEN} ${EXTENDED_SOURCE}") set_property(SOURCE rs_hinterp.F APPEND_STRING PROPERTY COMPILE_FLAGS "${BYTERECLEN} ${EXTENDED_SOURCE}") set_property(SOURCE time_ave.F APPEND_STRING PROPERTY COMPILE_FLAGS "${BYTERECLEN} ${EXTENDED_SOURCE}") +set_property(SOURCE Create_TEM_Diag.F APPEND_STRING PROPERTY COMPILE_FLAGS "${BYTERECLEN} ${EXTENDED_SOURCE}") set_property(SOURCE convert_aerosols.F APPEND_STRING PROPERTY COMPILE_FLAGS "${FREAL8} ${BYTERECLEN} ${EXTENDED_SOURCE}") set_property(SOURCE PWSSSP.F APPEND_STRING PROPERTY COMPILE_FLAGS "${FREAL8} ${EXTENDED_SOURCE} ${BYTERECLEN}") set_property(SOURCE alias.F APPEND_STRING PROPERTY COMPILE_FLAGS "${BYTERECLEN}") @@ -91,6 +94,7 @@ target_link_libraries (flat2hdf.x post_nompi GMAO_gfio_r4) target_link_libraries (gg2fv.x post) target_link_libraries (rs2hdf.x post) target_link_libraries (time_ave.x ${this} GMAO_gfio_r4 MAPL ) +target_link_libraries (Create_TEM_Diag.x ${this} GMAO_gfio_r4 ) target_link_libraries (3CH_mpi.x ${this} GMAO_gfio_r4 MPI::MPI_Fortran) target_link_libraries (rsg3_vinterp.x ${this} GMAO_hermes) diff --git a/post/Create_TEM_Diag b/post/Create_TEM_Diag new file mode 100755 index 00000000..31e3064a --- /dev/null +++ b/post/Create_TEM_Diag @@ -0,0 +1,102 @@ +#!/bin/csh -f + +if ( ! $?GEOSUTIL ) then + echo " " + echo Environment variable GEOSUTIL must be defined before use! + echo Set GEOSUTIL to the location of the \'install\' directory + echo " " + exit 1 +endif + +if( $#argv == 0 ) then + $GEOSUTIL/bin/Create_TEM_Diag.x + exit +endif + +if( "$argv[1]" == "-h" ) then + $GEOSUTIL/bin/Create_TEM_Diag.x + exit +else + $GEOSUTIL/bin/Create_TEM_Diag.x $argv + + if( -e basename.output ) /bin/rm basename.output + + # Create DDF File for TEM_Diags + # ----------------------------- + set TEM_Diags = `/bin/ls -1 *.TEM_Diag.monthly.* | grep -v clim | grep -v partial | grep -v ddf` + set num = $#TEM_Diags + + echo " " + echo "TEM_Diag files: $TEM_Diags" + echo "Number: $num" + echo " " + + if( $num != 0 ) then + + # Extract Extension (last node) from Monthly Mean File + # ---------------------------------------------------- + @ m = 1 + while( $m != 0 ) + set dum = `echo $TEM_Diags[1] | cut -d "." -f$m` + if( .$dum != . ) then + set ext = $dum + @ m = $m + 1 + else + @ m = 0 + endif + end + echo "Extension: $ext" + + # Extract EXPID from Monthly Mean File + # ------------------------------------ + set expid = `echo $TEM_Diags[1] | cut -d "." -f1` + @ m = 2 + while( $m != 0 ) + set dum = `echo $TEM_Diags[1] | cut -d "." -f$m` + if( .$dum != .TEM_Diag ) then + set expid = `echo ${expidr}.$dum` + @ m = $m + 1 + else + @ date_node = $m + 2 + @ m = 0 + endif + end + echo "EXPID: $expid" + echo " " + + set edate = `echo $TEM_Diags[$num] | cut -d "." -f${date_node}` + set eyear = `echo $edate | cut -c1-4` + set emonth = `echo $edate | cut -c5-6` + + set date = `echo $TEM_Diags[1] | cut -d "." -f${date_node}` + set year = `echo $date | cut -c1-4` + set month = `echo $date | cut -c5-6` + + set expdsc = "TEM Diagnostics" + set timinc = "mo" + + if( $month == "01" ) set MON = JAN + if( $month == "02" ) set MON = FEB + if( $month == "03" ) set MON = MAR + if( $month == "04" ) set MON = APR + if( $month == "05" ) set MON = MAY + if( $month == "06" ) set MON = JUN + if( $month == "07" ) set MON = JUL + if( $month == "08" ) set MON = AUG + if( $month == "09" ) set MON = SEP + if( $month == "10" ) set MON = OCT + if( $month == "11" ) set MON = NOV + if( $month == "12" ) set MON = DEC + + @ nmonths = ( $eyear - $year ) * 12 + ( $emonth - $month ) + 1 + + /bin/rm -f ${expid}.TEM_Diag.monthly.ddf + echo DSET ^$expid.TEM_Diag.monthly.%y4%m2.$ext > ${expid}.TEM_Diag.monthly.ddf + echo TITLE $expdsc >> ${expid}.TEM_Diag.monthly.ddf + echo OPTIONS template >> ${expid}.TEM_Diag.monthly.ddf + echo TDEF time $nmonths LINEAR 00:00Z01$MON$year 1$timinc >> ${expid}.TEM_Diag.monthly.ddf + + endif + +endif + diff --git a/post/Create_TEM_Diag.F b/post/Create_TEM_Diag.F new file mode 100644 index 00000000..d354914e --- /dev/null +++ b/post/Create_TEM_Diag.F @@ -0,0 +1,212 @@ + program main + implicit none + +! ********************************************************************** +! ********************************************************************** +! **** **** +! **** Program to create TEM_Diag data from Monthly Mean Prog File **** +! **** **** +! ********************************************************************** +! ********************************************************************** + + integer im,jm,lm,nt + integer nymd,nhms + + integer mem_unit,loc,ios + logical lopen + character*256 cbuff1 + character*1 cbuff2(256) + equivalence ( cbuff1,cbuff2 ) + character*256 expid + character*1 char + + character*256 filename + integer id,nvars,nsecf + integer ntime,ngatts,rc,timinc + real undef + integer i, m, n, TEM_NVARS, nfiles + + character*256 title + character*256 source + character*256 contact + character*256 levunits + character*256, allocatable :: vname(:) + character*256, allocatable :: vtitle(:) + character*256, allocatable :: vunits(:) + + real, allocatable :: TEM_VARS(:,:,:,:) + real, allocatable :: lat(:) + real, allocatable :: lon(:) + real, allocatable :: lev(:) + real, allocatable :: vrange(:,:) + real, allocatable :: prange(:,:) + integer, allocatable :: yymmdd(:) + integer, allocatable :: hhmmss(:) + integer, allocatable :: kmvar(:) + + integer nargs, iargc + character*256, allocatable :: arg(:) + character*256, allocatable :: fname(:) + +! ********************************************************************** + + nargs = iargc() + if( nargs.eq.0 ) then + call usage() + else + allocate ( arg(nargs) ) + do n=1,nargs + call getarg(n,arg(n)) + enddo + endif + +! ********************************************************************** +! **** Store INPUT into variables **** +! ********************************************************************** + + expid = 'NULL' + + do n=1,nargs + call getarg(n,arg(n)) + if( trim(arg(n)).eq.'-expid' ) expid = arg(n+1) + + if( trim(arg(n)).eq.'-files' ) then + nfiles = 1 + read(arg(n+nfiles),fmt='(a1)') char + do while (char.ne.'-' .and. n+nfiles.ne.nargs ) + nfiles = nfiles+1 + read(arg(n+nfiles),fmt='(a1)') char + enddo + if( char.eq.'-' ) nfiles = nfiles-1 + allocate ( fname(nfiles) ) + do m=1,nfiles + fname(m) = arg(n+m) + enddo + endif + + enddo + +! ********************************************************************** +! **** Extract EXPID from Filename if needed **** +! ********************************************************************** + + filename = fname(1) + + if( trim(expid).eq.'NULL' ) then + cbuff1 = '/usr/bin/basename ' // trim(filename) // ' > basename.output' + call execute_command_line ( trim(cbuff1) ) + + i = 1 + mem_unit = 0 + do while( mem_unit == 0 ) + if( i /= 5 .and. i /= 6 .and. i <= 99 ) then + inquire( unit = i, opened = lopen, iostat = ios ) + if( ios == 0 ) then + if( .not. lopen ) mem_unit = i + endif + endif + i = i+1 + enddo + + open (mem_unit,file='basename.output',form='formatted',access='sequential') + read (mem_unit,'(a)') cbuff1 + close(mem_unit) + + expid = '' + loc = 1 + char = cbuff2(loc) + do while (char.ne.'.') + expid = trim(expid) // char + loc = loc+1 + char = cbuff2(loc) + enddo + endif + +! ********************************************************************** +! **** Read Monthly Mean Prog File(s) **** +! ********************************************************************** + + do n=1,nfiles + + filename = fname(n) + + call gfio_open ( trim(filename),1,id,rc ) + call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) + + TEM_NVARS = 7 + + allocate ( lon(im) ) + allocate ( lat(jm) ) + allocate ( lev(lm) ) + allocate ( yymmdd(ntime) ) + allocate ( hhmmss(ntime) ) + allocate ( vname(nvars) ) + allocate ( vtitle(nvars) ) + allocate ( vunits(nvars) ) + allocate ( kmvar(nvars) ) + allocate ( vrange(2,nvars) ) + allocate ( prange(2,nvars) ) + + call gfio_inquire ( id,im,jm,lm,ntime,nvars, + . title,source,contact,undef, + . lon,lat,lev,levunits, + . yymmdd,hhmmss,timinc, + . vname,vtitle,vunits,kmvar, + . vrange,prange,rc ) + + nymd = yymmdd(1) + nhms = hhmmss(1) + + allocate ( TEM_VARS(im,jm,lm,TEM_NVARS) ) + + call gfio_getvar ( id,'U' ,nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,1),rc ) + call gfio_getvar ( id,'V' ,nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,2),rc ) + call gfio_getvar ( id,'T' ,nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,3),rc ) + call gfio_getvar ( id,'OMEGA',nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,4),rc ) + call gfio_getvar ( id,'USVS' ,nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,5),rc ) + call gfio_getvar ( id,'USWS' ,nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,6),rc ) + call gfio_getvar ( id,'VSTS' ,nymd,nhms,im,jm,1,lm,TEM_VARS(1,1,1,7),rc ) + + call gfio_close ( id,rc ) + +! ********************************************************************** +! **** Create TEM_Diag File **** +! ********************************************************************** + + call TEM( TEM_VARS,im,jm,lm,TEM_NVARS,expid,lon,lat,lev,levunits,nymd,nhms,timinc,undef ) + + deallocate ( lon ) + deallocate ( lat ) + deallocate ( lev ) + deallocate ( yymmdd ) + deallocate ( hhmmss ) + deallocate ( vname ) + deallocate ( vtitle ) + deallocate ( vunits ) + deallocate ( kmvar ) + deallocate ( vrange ) + deallocate ( prange ) + deallocate ( TEM_VARS ) + + enddo + +! ********************************************************************** + + stop + end + + subroutine usage() + write(6,100) + 100 format( "Usage: " ,/ + . ,/ + . " Create_TEM_Diag.x -files filename(s)" ,/ + . " [-expid expid]" ,/ + . ,/ + . "where:" ,/ + . ,/ + . " filename(s): Filename(s) containing pressure level data for U,V,T,OMEGA,USVS,USWS,VSTS" ,/ + . " expid: Optional EXPID to use for output filename (Default: 1st node of input filename[1]" ,/ + . ) + call exit(7) + end subroutine usage + diff --git a/post/TEM.F90 b/post/TEM.F90 new file mode 100644 index 00000000..12049ba1 --- /dev/null +++ b/post/TEM.F90 @@ -0,0 +1,1575 @@ + subroutine TEM( TEM_VARS,im,jm,lm,TEM_NVARS, & + EXPID,lon,lat,lev,levunits,nymd,nhms,timinc,undef ) + implicit none + character*6 date + + integer id,nvars,precision,rc + integer im,jm,lm,TEM_NVARS + integer nymd,nhms,ntime,timinc + integer i,j,L,imz + real undef + + real lat(jm) + real lon(im) + real lev(lm) + + real TEM_VARS(im,jm,lm,TEM_NVARS) + + character*256 filename + character*256 title + character*256 EXPID + character*256 output + character*256 source + character*256 contact + character*256 levunits + character*256, allocatable :: vname(:) + character*256, allocatable :: vtitle(:) + character*256, allocatable :: vunits(:) + + real, allocatable :: vrange(:,:) + real, allocatable :: prange(:,:) + integer, allocatable :: lmvar(:) + + real, allocatable :: pl3d(:,:,:) + real, allocatable :: pk3d(:,:,:) + real, allocatable :: u3d(:,:,:) + real, allocatable :: v3d(:,:,:) + real, allocatable :: t3d(:,:,:) + real, allocatable :: th3d(:,:,:) + real, allocatable :: omega3d(:,:,:) + real, allocatable :: upvp3d(:,:,:) + real, allocatable :: upwp3d(:,:,:) + real, allocatable :: vptp3d(:,:,:) + real, allocatable :: vpthp3d(:,:,:) + + real, allocatable :: uz(:,:) + real, allocatable :: vz(:,:) + real, allocatable :: tz(:,:) + real, allocatable :: wz(:,:) + real, allocatable :: omg(:,:) + real, allocatable :: thz(:,:) + real, allocatable :: upvpz(:,:) + real, allocatable :: upwpz(:,:) + real, allocatable :: vptpz(:,:) + real, allocatable :: vpthpz(:,:) + real, allocatable :: pl(:,:) + real, allocatable :: pk(:,:) + real, allocatable :: strm(:,:) + real, allocatable :: res(:,:) + real, allocatable :: vstar(:,:) + real, allocatable :: veddy(:,:) + real, allocatable :: wstar(:,:) + real, allocatable :: wmean(:,:) + real, allocatable :: weddy(:,:) + real, allocatable :: psi1(:,:) ! Residual Mass StreamFunction (Method 1) + real, allocatable :: psi2(:,:) ! Residual Mass StreamFunction (Method 2) + real, allocatable :: psim(:,:) ! Mean Mass StreamFunction + real, allocatable :: epfy(:,:) ! Eliassen-Palm Flux in Northward Direction + real, allocatable :: epfz(:,:) ! Eliassen-Palm Flux in Upward Direction + real, allocatable :: epfdiv(:,:) ! Eliassen-Palm Flux Divergence + real, allocatable :: vstr(:,:) + real, allocatable :: vstar2(:,:) + real, allocatable :: wstar2(:,:) + real, allocatable :: wmean2(:,:) + real, allocatable :: weddy2(:,:) + real*4, allocatable :: dum(:) + + real,allocatable :: upvp (:,:) + real,allocatable :: upwp (:,:) + real,allocatable :: dudp (:,:) + real,allocatable :: dudphi(:,:) + real,allocatable :: psie (:,:) + real,allocatable :: dfdphi(:,:) + real,allocatable :: dfdp (:,:) + real,allocatable :: plz (:,:) + real,allocatable :: delp (:,:) + + integer n,nt + +! Read data from grads.fwrite: Data is written BOT (L=1) to TOP (L=LM) +! -------------------------------------------------------------------- + + ntime = 1 + + print * + print *, 'Creating TEM_Diag Data for EXPID: ',trim(EXPID) + print *, '--------------------------------- ' + print * + print *, 'IM: ',im + print *, 'JM: ',jm + print *, 'LM: ',lm + print * + print *, ' NYMD: ',nymd + print *, ' NHMS: ',nhms + print * + + allocate ( pl3d(im,jm,lm) ) + allocate ( pk3d(im,jm,lm) ) + + allocate ( u3d(im,jm,lm) ) + allocate ( v3d(im,jm,lm) ) + allocate ( t3d(im,jm,lm) ) + allocate ( th3d(im,jm,lm) ) + allocate ( upvp3d(im,jm,lm) ) + allocate ( upwp3d(im,jm,lm) ) + allocate ( vptp3d(im,jm,lm) ) + allocate ( vpthp3d(im,jm,lm) ) + allocate ( omega3d(im,jm,lm) ) + + ! Initialize Variables from BOT (L=1) to TOP (L=LM) + ! ------------------------------------------------- + do L=1,lm + u3d(:,:,L) = TEM_VARS(:,:,L,1) + v3d(:,:,L) = TEM_VARS(:,:,L,2) + t3d(:,:,L) = TEM_VARS(:,:,L,3) + omega3d(:,:,L) = TEM_VARS(:,:,L,4) + upvp3d(:,:,L) = TEM_VARS(:,:,L,5) + upwp3d(:,:,L) = TEM_VARS(:,:,L,6) + vptp3d(:,:,L) = TEM_VARS(:,:,L,7) + enddo + + ! Define Pressure Variables from BOT (L=1) to TOP (L=LM) + ! ------------------------------------------------------ + do L=1,lm + do j=1,jm + do i=1,im + pl3d(i,j,L) = lev(L) ! mb + pk3d(i,j,L) = pl3d(i,j,L)**(2.0/7.0) ! mb**kappa + enddo + enddo + enddo + + do L=1,lm + ! print *, 'L: ',L,' LEV: ',pl3d(1,jm/2,L),' T: ',t3d(1,jm/2,L) + do j=1,jm + do i=1,im + if( t3d(i,j,L).ne.undef ) then + th3d(i,j,L) = t3d(i,j,L)/pk3d(i,j,L) ! K/(mb**kappa) + else + th3d(i,j,L) = undef + endif + if( vptp3d(i,j,L).ne.undef ) then + vpthp3d(i,j,L) = vptp3d(i,j,L)/pk3d(i,j,L) ! m/sec K / (mp**kappa) + else + vpthp3d(i,j,L) = undef + endif + enddo + enddo + enddo + + allocate( uz(jm,lm) ) + allocate( vz(jm,lm) ) + allocate( tz(jm,lm) ) + allocate( wz(jm,lm) ) + allocate( omg(jm,lm) ) + allocate( thz(jm,lm) ) + allocate( upvpz(jm,lm) ) + allocate( upwpz(jm,lm) ) + allocate( vptpz(jm,lm) ) + allocate( vpthpz(jm,lm) ) + allocate( pl(jm,lm) ) + allocate( pk(jm,lm) ) + allocate( strm(jm,lm) ) + allocate( res(jm,lm) ) + allocate( vstar(jm,lm) ) + allocate( veddy(jm,lm) ) + allocate( wstar(jm,lm) ) + allocate( wmean(jm,lm) ) + allocate( weddy(jm,lm) ) + allocate( psi1(jm,lm) ) + allocate( psi2(jm,lm) ) + allocate( psim(jm,lm) ) + allocate( epfy(jm,lm) ) + allocate( epfz(jm,lm) ) + allocate( epfdiv(jm,lm) ) + allocate( vstr(jm,lm) ) + allocate( dum(jm) ) + allocate( vstar2(jm,lm) ) + allocate( wstar2(jm,lm) ) + allocate( wmean2(jm,lm) ) + allocate( weddy2(jm,lm) ) + + allocate( upvp (jm,LM) ) + allocate( upwp (jm,LM) ) + allocate( dudp (jm,LM) ) + allocate( dudphi(jm,LM) ) + allocate( psie (jm,LM) ) + allocate( dfdphi(jm,LM) ) + allocate( dfdp (jm,LM) ) + allocate( plz (jm,LM) ) + allocate( delp (jm,LM) ) + + call zonal ( v3d,im,jm,lm, vz,undef) + call zonal ( t3d,im,jm,lm, tz,undef) + call zonal ( vptp3d,im,jm,lm, vptpz,undef) + + call zonal ( pl3d,im,jm,lm, pl,undef) + call zonal ( pk3d,im,jm,lm, pk,undef) + call zonal ( u3d,im,jm,lm, uz,undef) + call zonal ( upvp3d,im,jm,lm, upvpz,undef) + call zonal ( upwp3d,im,jm,lm, upwpz,undef) + call zonal (omega3d,im,jm,lm, omg,undef) + + call zonal ( th3d,im,jm,lm, thz,undef) + call zonal (vpthp3d,im,jm,lm,vpthpz,undef) + +! ------------------------------------------------------------------------------ + + do L=1,lm + do j=1,jm + if( abs(tz(j,L)-undef).gt.0.1 ) then + thz(j,L) = tz(j,L)/pk(j,L) ! K/mb**kappa + else + thz(j,L) = undef + endif + if( abs(vptpz(j,L)-undef).gt.0.1 ) then + vpthpz(j,L) = vptpz(j,L)/pk(j,L) ! m/sec K/mb**kappa + else + vpthpz(j,L) = undef + endif + enddo + enddo + +! Compute Meridional Streamfunction +! --------------------------------- + call stream ( vz,pl,jm,lm,strm,undef ) + + call make_psi ( uz,vz,thz,omg,upvpz,upwpz,vpthpz,pl,jm,lm,psi1,psi2,psim,epfy,epfz,epfdiv, & + undef,upvp,upwp,dudp,dudphi,psie,dfdphi,dfdp,vstar2,veddy,wstar2,wmean2,weddy2,plz,delp) + +! Compute Mean Vertical Velocity from Continuity +! ---------------------------------------------- + call make_w ( vz,pl,jm,lm,wz,undef ) + +! Compute Residual Circulation using wmean2 from model data rather than wz from continuity +! ---------------------------------------------------------------------------------------- + call residual ( vz,vpthpz,thz,wmean2,pl,jm,lm,res,vstar,wstar,wmean,weddy,undef ) + + +! ********************************************************************** +! **** Write TEM_Diag Monthly Output File **** +! ********************************************************************** + + title = 'Transformed Eulerian Mean (TEM) Diagnostics' + source = 'Goddard Modeling and Assimilation Office, NASA/GSFC' + contact = 'data@gmao.gsfc.nasa.gov' + + imz = 1 + nvars = 26 + + allocate ( vname(nvars) ) + allocate ( vtitle(nvars) ) + allocate ( vunits(nvars) ) + allocate ( lmvar(nvars) ) + allocate ( vrange(2,nvars) ) + allocate ( prange(2,nvars) ) + + vrange(:,:) = undef + prange(:,:) = undef + lmvar (:) = lm + vunits(:) = 'unknown' + + vname(01) = 'str' ; vtitle(01) = 'Streamfunction' + vname(02) = 'res' ; vtitle(02) = 'Residual Circulation' + vname(03) = 'vstar' ; vtitle(03) = 'vstar' + vname(04) = 'wstar' ; vtitle(04) = 'wstar' + vname(05) = 'wmean' ; vtitle(05) = 'wmean' + vname(06) = 'weddy' ; vtitle(06) = 'weddy' + vname(07) = 'psi1' ; vtitle(07) = 'Res1 Streamfunction' + vname(08) = 'psi2' ; vtitle(08) = 'Res2 Streamfunction' + vname(09) = 'psim' ; vtitle(09) = 'Mass Streamfunction' + vname(10) = 'epfy' ; vtitle(10) = 'Eliassen-Palm flux y' + vname(11) = 'epfz' ; vtitle(11) = 'Eliassen-Palm flux z' + vname(12) = 'epfdiv' ; vtitle(12) = 'Eliassen-Palm flux Divergence' + vname(13) = 'upvp' ; vtitle(13) = 'Uprime Vprime' + vname(14) = 'upwp' ; vtitle(14) = 'Uprime Omegaprime' + vname(15) = 'dudp' ; vtitle(15) = 'DuDp' + vname(16) = 'dudphi' ; vtitle(16) = 'DuDphi' + vname(17) = 'psie' ; vtitle(17) = 'Eddy Streamfunction' + vname(18) = 'dfdphi' ; vtitle(18) = 'DfDphi' + vname(19) = 'dfdp' ; vtitle(19) = 'DfDp' + vname(20) = 'plz' ; vtitle(20) = 'Pressure' + vname(21) = 'delp' ; vtitle(21) = 'Pressure Thickness' + vname(22) = 'vstar2' ; vtitle(22) = 'Alternate vstar Calculation' + vname(23) = 'veddy2' ; vtitle(23) = 'Alternate veddy Calculation' + vname(24) = 'wstar2' ; vtitle(24) = 'Alternate wstar Calculation' + vname(25) = 'wmean2' ; vtitle(25) = 'Alternate wmean Calculation' + vname(26) = 'weddy2' ; vtitle(26) = 'Alternate weddy Calculation' + + precision = 1 ! 64-bit + precision = 0 ! 32-bit + + write(date,1000) nymd/100 + 1000 format(i6.6) + + filename = trim(EXPID) // '.TEM_Diag.monthly.' // date // '.nc4' + + call GFIO_Create ( trim(filename), trim(title), source, contact, undef, & + imz, jm, lm, lon(1), lat, lev, levunits, & + nymd, nhms, timinc, & + nvars, vname, vtitle, vunits, lmvar, & + vrange, prange, precision, & + id, rc ) + + call Gfio_putVar ( id,trim(vname(01)),nymd,nhms,imz,jm,1,lm,strm ,rc ) ; print *, 'Writing variable: ',trim(vname(01)) + call Gfio_putVar ( id,trim(vname(02)),nymd,nhms,imz,jm,1,lm,res ,rc ) ; print *, 'Writing variable: ',trim(vname(02)) + call Gfio_putVar ( id,trim(vname(03)),nymd,nhms,imz,jm,1,lm,vstar ,rc ) ; print *, 'Writing variable: ',trim(vname(03)) + call Gfio_putVar ( id,trim(vname(04)),nymd,nhms,imz,jm,1,lm,wstar ,rc ) ; print *, 'Writing variable: ',trim(vname(04)) + call Gfio_putVar ( id,trim(vname(05)),nymd,nhms,imz,jm,1,lm,wmean ,rc ) ; print *, 'Writing variable: ',trim(vname(05)) + call Gfio_putVar ( id,trim(vname(06)),nymd,nhms,imz,jm,1,lm,weddy ,rc ) ; print *, 'Writing variable: ',trim(vname(06)) + call Gfio_putVar ( id,trim(vname(07)),nymd,nhms,imz,jm,1,lm,psi1 ,rc ) ; print *, 'Writing variable: ',trim(vname(07)) + call Gfio_putVar ( id,trim(vname(08)),nymd,nhms,imz,jm,1,lm,psi2 ,rc ) ; print *, 'Writing variable: ',trim(vname(08)) + call Gfio_putVar ( id,trim(vname(09)),nymd,nhms,imz,jm,1,lm,psim ,rc ) ; print *, 'Writing variable: ',trim(vname(09)) + call Gfio_putVar ( id,trim(vname(10)),nymd,nhms,imz,jm,1,lm,epfy ,rc ) ; print *, 'Writing variable: ',trim(vname(10)) + call Gfio_putVar ( id,trim(vname(11)),nymd,nhms,imz,jm,1,lm,epfz ,rc ) ; print *, 'Writing variable: ',trim(vname(11)) + call Gfio_putVar ( id,trim(vname(12)),nymd,nhms,imz,jm,1,lm,epfdiv ,rc ) ; print *, 'Writing variable: ',trim(vname(12)) + call Gfio_putVar ( id,trim(vname(13)),nymd,nhms,imz,jm,1,lm,upvp ,rc ) ; print *, 'Writing variable: ',trim(vname(13)) + call Gfio_putVar ( id,trim(vname(14)),nymd,nhms,imz,jm,1,lm,upwp ,rc ) ; print *, 'Writing variable: ',trim(vname(14)) + call Gfio_putVar ( id,trim(vname(15)),nymd,nhms,imz,jm,1,lm,dudp ,rc ) ; print *, 'Writing variable: ',trim(vname(15)) + call Gfio_putVar ( id,trim(vname(16)),nymd,nhms,imz,jm,1,lm,dudphi ,rc ) ; print *, 'Writing variable: ',trim(vname(16)) + call Gfio_putVar ( id,trim(vname(17)),nymd,nhms,imz,jm,1,lm,psie ,rc ) ; print *, 'Writing variable: ',trim(vname(17)) + call Gfio_putVar ( id,trim(vname(18)),nymd,nhms,imz,jm,1,lm,dfdphi ,rc ) ; print *, 'Writing variable: ',trim(vname(18)) + call Gfio_putVar ( id,trim(vname(19)),nymd,nhms,imz,jm,1,lm,dfdp ,rc ) ; print *, 'Writing variable: ',trim(vname(19)) + call Gfio_putVar ( id,trim(vname(20)),nymd,nhms,imz,jm,1,lm,plz ,rc ) ; print *, 'Writing variable: ',trim(vname(20)) + call Gfio_putVar ( id,trim(vname(21)),nymd,nhms,imz,jm,1,lm,delp ,rc ) ; print *, 'Writing variable: ',trim(vname(21)) + call Gfio_putVar ( id,trim(vname(22)),nymd,nhms,imz,jm,1,lm,vstar2 ,rc ) ; print *, 'Writing variable: ',trim(vname(22)) + call Gfio_putVar ( id,trim(vname(23)),nymd,nhms,imz,jm,1,lm,veddy ,rc ) ; print *, 'Writing variable: ',trim(vname(23)) + call Gfio_putVar ( id,trim(vname(24)),nymd,nhms,imz,jm,1,lm,wstar2 ,rc ) ; print *, 'Writing variable: ',trim(vname(24)) + call Gfio_putVar ( id,trim(vname(25)),nymd,nhms,imz,jm,1,lm,wmean2 ,rc ) ; print *, 'Writing variable: ',trim(vname(25)) + call Gfio_putVar ( id,trim(vname(26)),nymd,nhms,imz,jm,1,lm,weddy2 ,rc ) ; print *, 'Writing variable: ',trim(vname(26)) + print * + + call gfio_close ( id,rc ) + +#if 0 +! Write Grads Binary Data +! ----------------------- + call GLAWRT (strm ,jm,LM,20,undef) + call GLAWRT (res ,jm,LM,20,undef) + call GLAWRT (vstar ,jm,LM,20,undef) + call GLAWRT (wstar ,jm,LM,20,undef) + call GLAWRT (wmean ,jm,LM,20,undef) + call GLAWRT (weddy ,jm,LM,20,undef) + call GLAWRT (psi1 ,jm,LM,20,undef) + call GLAWRT (psi2 ,jm,LM,20,undef) + call GLAWRT (psim ,jm,LM,20,undef) + call GLAWRT (epfy ,jm,LM,20,undef) + call GLAWRT (epfz ,jm,LM,20,undef) + call GLAWRT (epfdiv,jm,LM,20,undef) + + call GLAWRT (upvp ,jm,LM,20,undef) + call GLAWRT (upwp ,jm,LM,20,undef) + call GLAWRT (dudp ,jm,LM,20,undef) + call GLAWRT (dudphi,jm,LM,20,undef) + call GLAWRT (psie ,jm,LM,20,undef) + call GLAWRT (dfdphi,jm,LM,20,undef) + call GLAWRT (dfdp ,jm,LM,20,undef) + call GLAWRT (plz ,jm,LM,20,undef) + call GLAWRT (delp ,jm,LM,20,undef) + call GLAWRT (vstar2,jm,LM,20,undef) + call GLAWRT (veddy ,jm,LM,20,undef) + call GLAWRT (wstar2,jm,LM,20,undef) + call GLAWRT (wmean2,jm,LM,20,undef) + call GLAWRT (weddy2,jm,LM,20,undef) + +! Write Grads Control File +! ------------------------ + output = '^residual.' // 'test' // '.data' + output = '^fort.20' + + write(30,101) trim(output),trim(title),undef,jm,lat(1),2.0*abs(lat(1))/(jm-1),lm,lev(1:lm) + do L=1,lm + print *, 'Pressure = ',lev(L) + enddo + print * + write(30,103) '00Z01DEC2012',lm,lm,lm,lm,lm,lm,lm,lm,lm,lm,lm,lm, & + lm,lm,lm,lm,lm,lm,lm,lm,lm,lm,lm,lm,lm,lm + 101 format('dset ',a,/, & + 'title ',a,/, & + 'options sequential big_endian',/, & + 'undef ',e15.6,/, & + 'xdef 1 linear -180 1',/, & + 'ydef ',i4,' linear ',f8.3,2x,f8.3,/, & + 'zdef ',i3,' levels ',999(f8.3,1x)) + 102 format(10x,f8.3) + 103 format('tdef ','1 linear ',a,' 1mo',/, & + 'vars 26',/, & + 'str ',i3,' 0 Streamfunction',/, & + 'res ',i3,' 0 Residual Circulation',/, & + 'vstar ',i3,' 0 Vstar',/, & + 'wstar ',i3,' 0 wstar',/, & + 'wmean ',i3,' 0 wmean ',/, & + 'weddy ',i3,' 0 weddy ',/, & + 'psi1 ',i3,' 0 Res1 streamfunction ',/, & + 'psi2 ',i3,' 0 Res2 streamfunction ',/, & + 'psim ',i3,' 0 Mass streamfunction ',/, & + 'epfy ',i3,' 0 Eliassen-Palm flux y',/, & + 'epfz ',i3,' 0 Eliassen-Palm flux z',/, & + 'epfdiv ',i3,' 0 Eliassen-Palm flux Divergence',/, & + 'upvp ',i3,' 0 Uprime Vprim',/, & + 'upwp ',i3,' 0 Uprime Omegaprime',/, & + 'dudp ',i3,' 0 DuDp',/, & + 'dudphi ',i3,' 0 DuDphi',/, & + 'psie ',i3,' 0 Eddy Streamfunction',/, & + 'dfdphi ',i3,' 0 DfDphi',/, & + 'dfdp ',i3,' 0 DfDp',/, & + 'plz ',i3,' 0 Pressure',/, & + 'delp ',i3,' 0 Pressure Thickness',/, & + 'vstar2 ',i3,' 0 vstar2',/, & + 'veddy2 ',i3,' 0 veddy2 ',/, & + 'wstar2 ',i3,' 0 wstar2',/, & + 'wmean2 ',i3,' 0 wmean2',/, & + 'weddy2 ',i3,' 0 weddy2',/, & + 'endvars') +#endif + + return + end + + ! ************************************************************************************************************ + + subroutine zonal (q,im,jm,lm,qz,undef) + implicit none + integer im,jm,lm + integer i,j,L,isum + real q(im,jm,lm), qz(jm,lm) + real undef + real*8 qsum + logical defined + + do L=1,lm + do j=1,jm + qsum = 0.0D0 + isum = 0 + do i=1,im + if( defined(q(i,j,L),undef) ) then + qsum = qsum + q(i,j,L) + isum = isum + 1 + endif + enddo + if( isum.ne.0 ) then + qz(j,L) = qsum/isum + else + qz(j,L) = undef + endif + enddo + enddo + + return + end + + ! ************************************************************************************************************ + + SUBROUTINE GLAWRT (A, JM,LM, KTP, undef) + real A (JM,LM) + real*4 TEM (JM), undef + logical defined + DO L=1,LM + DO J=1,JM + if( defined(a(J,L),undef) ) then + if( abs(a(J,L)).gt.1.e-20 ) then + TEM (J) = A(J,L) + else + TEM (J) = 0. + endif + else + TEM (J) = undef + endif + ENDDO + WRITE(KTP) TEM + ENDDO + RETURN + END + + subroutine make_psi( u0,v0,th0,w0,upvp0,upwp0,vpthp0,p0,jm,lm,psi1,psi2,psim,epfy,epfz,epfdiv,undef, & + upvp,upwp,dudp,dudphi,psie,dfdphi,dfdp,vstar,veddy,wstar,wmean,weddy,p,delp) + use MAPL_ConstantsMod + implicit none + integer j,k,L,jm,lm + real undef,dphi,a,g,pi,phi,pk0,H + logical defined + + real th0(jm,lm), th(jm,lm) + real upvp0(jm,lm), upvp(jm,lm) + real upwp0(jm,lm), upwp(jm,lm) + real vpthp0(jm,lm), vpthp(jm,lm) + real u0(jm,lm), u(jm,lm) + real v0(jm,lm), v(jm,lm) + real p0(jm,lm), p(jm,lm) + real w0(jm,lm), w(jm,lm) + + real psi1(jm,lm) + real psi2(jm,lm) + real psim(jm,lm) + real epfy(jm,lm) + real epfz(jm,lm) + real epfdiv(jm,lm) + + real dudp(jm,lm) + real dfdp(jm,lm) + real dthdp(jm,lm) + real dudphi(jm,lm) + real dfdphi(jm,lm) + real psie(jm,lm) + real delp(jm,lm) + real veddy(jm,lm) + real vstar(jm,lm) + real weddy(jm,lm) + real wstar(jm,lm) + real wmean(jm,lm) + real stuff(jm,lm) + real dume(jm,0:lm) ! dummy edge + real the(jm,0:lm) ! theta_edge + real ple(jm,0:lm) ! p_edge + real pm(jm,1:lm) ! mid-point of p_edges + real ue(jm,0:lm) ! u_edge + real epfze(jm,0:lm) ! epfz_edge + real f(jm) + real dum(jm) + integer method + real airmw,runiv,cpd,rgas,akap, sum + + PARAMETER ( AIRMW = MAPL_AIRMW ) + PARAMETER ( RUNIV = MAPL_RUNIV ) + PARAMETER ( CPD = MAPL_CP ) + PARAMETER ( RGAS = RUNIV/AIRMW ) + PARAMETER ( AKAP = MAPL_KAPPA ) + +! Define Constants +! ---------------- + pi = 4.*atan(1.) + dphi = pi/(jm-1) + a = MAPL_RADIUS + g = MAPL_GRAV + H = 7000.0 + + Method = 0 + +! Invert level index (in order to be top=>bottom) +! ----------------------------------------------- + do L=1,lm + u(:,L) = u0(:,lm-L+1) ! m/sec + v(:,L) = v0(:,lm-L+1) ! m/sec + w(:,L) = w0(:,lm-L+1) ! Pa/sec + p(:,L) = p0(:,lm-L+1) ! mb + th(:,L) = th0(:,lm-L+1) ! K/mb**kappa + upvp(:,L) = upvp0(:,lm-L+1) ! m/sec m/sec + upwp(:,L) = upwp0(:,lm-L+1) ! m/sec Pa/sec + vpthp(:,L) = vpthp0(:,lm-L+1) ! m/sec K/mb**kappa + enddo + + pk0 = (1000.0)**(2.0/7.0) ! mb**kappa + + where( abs(p -undef).gt.0.1 ) ; p = p*100 ; endwhere + where( abs(th -undef).gt.0.1 ) ; th = th*pk0 ; endwhere + where( abs(vpthp-undef).gt.0.1 ) ; vpthp = vpthp*pk0 ; endwhere + +! Compute PLE Edge Values +! ----------------------- + ple(:,0) = max( 0.0, p(:,1) - 0.5*( p(:,2)-p(:,1) ) ) + do L=1,lm-1 + do j=1,jm + ple(j,L) = ( p(j,L+1)+ p(j,L) )*0.5 + enddo + enddo + ple(:,lm) = p(:,lm) + 0.5*( p(:,lm)-p(:,lm-1) ) + + do L=1,lm + delp(:,L) = ple(:,L)-ple(:,L-1) + enddo + +! Compute Mid-Point of PLE Edge Values +! ------------------------------------ + do L=1,lm + pm(:,L) = ( ple(:,L)+ple(:,L-1) )*0.5 + enddo + +! Compute Mass Streamfunction +! --------------------------- + pi = 4.*atan(1.) + dphi = pi/(jm-1) + a = MAPL_RADIUS + g = MAPL_GRAV + + do L=1,LM + dum(:) = 0.0 + do k=1,L + where( abs(v(:,k)-undef).gt.0.1 ) + dum(:) = dum(:) + v(:,k)*delp(:,k) + endwhere + enddo + do j=1,jm + phi = -pi/2 + (j-1)*dphi + psim(j,L) = 2*pi*a*cos(phi)/g * dum(j) + enddo + enddo + + sum = 0.0 + do k=1,lm + do j=1,jm + sum = sum + abs( psim(j,k) ) + enddo + enddo + if( sum.eq.0.0 ) psim = undef + +! Define Eddy Streamfunction: psie = vpthp/dthdp +! ----------------------------------------------- + + call map1_cubic( lm,p,th, lm+1,ple,the, jm, Method, undef) + call compute_d_dp( the,delp,jm,lm,undef,stuff ) + call map1_cubic( lm,pm,stuff, lm,p,dthdp, jm, Method, undef) + + do L=1,lm + do j=1,jm + if( defined(dthdp(j,L),undef) .and. & + defined(vpthp(j,L),undef) ) then + dthdp(j,L) = min( -0.003*pk0/100, dthdp(j,L) ) + psie(j,L) = vpthp(j,L) / dthdp(j,L) + else + psie(j,L) = undef + endif + enddo + enddo + +! Compute Veddy = D/Dp[ psie ] +! ---------------------------- + + call map1_cubic( lm,p,psie, lm+1,ple,dume, jm, Method, undef) + call compute_d_dp( dume,delp,jm,lm,undef,stuff ) + call map1_cubic( lm,pm,stuff, lm,p,veddy, jm, Method, undef) + +! Compute Vstar = v - veddy +! ------------------------- + do L=1,lm + do j=1,jm + if( defined( veddy(j,L),undef) .and. & + defined( v(j,L),undef) ) then + vstar(j,L) = v(j,L) - veddy(j,L) + else + vstar(j,L) = undef + endif + enddo + enddo + + +! Compute weddy = -d(psie*cos)/(a*cos*dphi) +! ----------------------------------------- + + do L=1,lm + do j=1,jm + if( defined( psie(j,L),undef ) ) then + stuff(j,L) = psie(j,L) + else + stuff(j,L) = undef + endif + enddo + enddo + + call compute_d_dphi( stuff,jm,lm,undef,weddy ) + + +! Compute wstar = w - weddy +! ------------------------- + do L=1,lm + do j=1,jm + if( defined( weddy(j,L),undef) .and. & + defined( w(j,L),undef) ) then + wstar(j,L) = w(j,L) + weddy(j,L) + else + wstar(j,L) = undef + endif + enddo + enddo + + do L=1,lm + do j=1,jm + if( defined( wstar(j,L),undef) ) then + wstar(j,L) = -H*wstar(j,L)/p(j,L) + endif + if( defined( w(j,L),undef) ) then + wmean(j,L) = -H*w(j,L)/p(j,L) + else + wmean(j,L) = undef + endif + if( defined( weddy(j,L),undef) ) then + weddy(j,L) = -H*weddy(j,L)/p(j,L) + endif + enddo + enddo + + +! Construct Residual Streamfunction from Vstar (Method 1) +! ------------------------------------------------------- + do L=1,LM + dum(:) = 0.0 + do k=1,L + where( abs(vstar(:,k)-undef).gt.0.1 ) + dum(:) = dum(:) + vstar(:,k)*delp(:,k) + endwhere + enddo + do j=1,jm + phi = -pi/2 + (j-1)*dphi + psi1(j,L) = 2*pi*a*cos(phi)/g * dum(j) + enddo + enddo + + sum = 0.0 + do k=1,lm + do j=1,jm + sum = sum + abs( psi1(j,k) ) + enddo + enddo + if( sum.eq.0.0 ) psi1 = undef + + +! Compute Residual Streamfunction (Method 2) +! ------------------------------------------ + do L=1,lm + do j=1,jm + phi = -pi/2 + (j-1)*dphi + if( defined(psie(j,L),undef) ) then + psi2(j,L) = 2*pi*a*cos(phi)/g * psie(j,L) + else + psi2(j,L) = undef + endif + enddo + enddo + + do L=1,lm + where( abs(psim(:,L)-undef).gt.0.1 .and. & + abs(psi2(:,L)-undef).gt.0.1 ) + psi2(:,L) = psim(:,L) - psi2(:,L) + elsewhere + psi2(:,L) = undef + endwhere + enddo + + +! Compute Eliassen-Palm Flux +! -------------------------- + do j=1,jm + phi = -pi/2 + (j-1)*dphi + f(j) = 2*MAPL_OMEGA*sin(phi) + enddo + + !------------------------- Compute du/dp -------------------------------- + + call map1_cubic( lm,p,u, lm+1,ple,ue, jm, Method, undef) + call compute_d_dp( ue,delp,jm,lm,undef,stuff ) + call map1_cubic( lm,pm,stuff, lm,p,dudp, jm, Method, undef) + + !--------------------- Compute d(u*cos)/(a*cos*dphi) --------------------- + + call compute_d_dphi( u,jm,lm,undef,dudphi ) + + !----------------------- Compute epfy & epfz ---------------------------- + + do L=1,lm + do j=1,jm + phi = -pi/2 + (j-1)*dphi + if( defined( dudp(j,L),undef) .and. & + defined( psie(j,L),undef) .and. & + defined( upvp(j,L),undef) ) then + epfy(j,L) = a*cos(phi)*( dudp(j,L)*psie(j,L) - upvp(j,L) ) + else + epfy(j,L) = undef + endif + if( defined( dudphi(j,L),undef) .and.& + defined( psie(j,L),undef) .and. & + defined( upwp(j,L),undef) ) then + epfz(j,L) = a*cos(phi)*( (f(j)-dudphi(j,L))*psie(j,L) - upwp(j,L) ) + else + epfz(j,L) = undef + endif + enddo + enddo + + !----------------------- Compute d(epfy*cos)/(a*cos*dphi) ----------------------- + + call compute_d_dphi( epfy,jm,lm,undef,dfdphi ) + + !------------------------- Compute d(epfz)/dp --------------------------- + + call map1_cubic( lm,p,epfz, lm+1,ple,epfze, jm, Method, undef) + call compute_d_dp( epfze,delp,jm,lm,undef,stuff ) + call map1_cubic( lm,pm,stuff, lm,p,dfdp, jm, Method, undef) + + !----------------------- Compute EPFlux Divergence ---------------------- + + do L=1,lm + do j=1,jm + if( defined( dfdp(j,L),undef) .and. & + defined( dfdphi(j,L),undef) ) then + epfdiv(j,L) = dfdphi(j,L) + dfdp(j,L) + else + epfdiv(j,L) = undef + endif + enddo + enddo + +! Invert Streamfunction for grads output (in order to be bottom=>top) +! ------------------------------------------------------------------- + + call flipz( psim ,jm,lm,1.0e-10 ,undef,'psim' ) + call flipz( psi1 ,jm,lm,1.0e-10*2.4892,undef,'psi1' ) + call flipz( psi2 ,jm,lm,1.0e-10*2.4892,undef,'psi2' ) + call flipz( epfy ,jm,lm,1.0 ,undef,'epfy' ) + call flipz( epfz ,jm,lm,1.0 ,undef,'epfz' ) + call flipz( epfdiv,jm,lm,1.0 ,undef,'epfdiv' ) + + call flipz( upvp ,jm,lm,1.0 ,undef,'upvp' ) + call flipz( upwp ,jm,lm,1.0 ,undef,'upwp' ) + call flipz( dudp ,jm,lm,1.0 ,undef,'dudp' ) + call flipz( dudphi,jm,lm,1.0 ,undef,'dudphi' ) + call flipz( psie ,jm,lm,1.0 ,undef,'psie' ) + call flipz( dfdphi,jm,lm,1.0 ,undef,'dfdphi' ) + call flipz( dfdp ,jm,lm,1.0 ,undef,'dfdp' ) + call flipz( p ,jm,lm,1.0 ,undef,'p' ) + call flipz( delp ,jm,lm,1.0 ,undef,'delp' ) + + call flipz( vstar ,jm,lm,1.0 ,undef,'vstar' ) + call flipz( veddy ,jm,lm,1.0 ,undef,'veddy' ) + call flipz( wstar ,jm,lm,1.0 ,undef,'wstar' ) + call flipz( wmean ,jm,lm,1.0 ,undef,'wmean' ) + call flipz( weddy ,jm,lm,1.0 ,undef,'weddy' ) + + return + end + + subroutine flipz( q,jm,lm,scale,undef,name ) + implicit none + character(*) name + integer j,L,jm,lm + real undef,scale + logical defined + real q(jm,lm) + real z(jm,lm) + do L=1,lm + do j=1,jm + if( isnan(q(j,LM-L+1)) ) then + print *, trim(name),' LM-L+1: ',LM-L+1,' J: ',j,' q: ',q(j,LM-L+1) + endif + if( defined(q(j,LM-L+1),undef) ) then + z(j,L) = q(j,LM-L+1)*scale + else + z(j,L) = undef + endif + enddo + enddo + do L=1,lm + q(:,L) = z(:,L) + enddo + + return + end + + subroutine stream ( v0,p0,jm,lm,s,undef ) + use MAPL_ConstantsMod + implicit none + integer j,k,L,jm,lm + real pi,dp,a,g,const,phi,undef + + real v(jm,lm), v0(jm,lm) + real s(jm,lm) + real p0(jm,lm), p(jm,lm) + real dum(jm), sum + + real ple(jm,0:lm) + real delp(jm, lm) + +! Invert VWND and P level index (in order to be top=>bottom) +! ---------------------------------------------------------- + do L=1,lm + p(:,L) = p0(:,lm-L+1) + v(:,L) = v0(:,lm-L+1) + enddo + +! Compute Edge Pressures and Thickness +! ------------------------------------ + ple(:,0) = max( 0.0, p(:,1) - 0.5*( p(:,2)-p(:,1) ) ) + do L=1,lm-1 + ple(:,L) = ( p(:,L)+ p(:,L+1) )*0.5 + enddo + ple(:,lm) = p(:,lm) + 0.5*( p(:,lm)-p(:,lm-1) ) + do L=1,lm + delp(:,L) = ple(:,L)-ple(:,L-1) + enddo + + pi = 4.*atan(1.) + dp = pi/(jm-1) + a = MAPL_RADIUS + g = MAPL_GRAV + + const = 2*pi*a/g * 1.0e-8 + + do k=1,lm + dum(:) = 0.0 + do L=1,k + do j=1,jm + phi = -pi/2+(j-1)*dp + if( abs(v(j,L)-undef).gt.0.1 ) then + dum(j) = dum(j) + v(j,L)*cos(phi)*delp(j,L) + endif + enddo + enddo + s(:,k) = dum(:)*const + enddo + + sum = 0.0 + do k=1,lm + do j=1,jm + sum = sum + abs( s(j,k) ) + enddo + enddo + if( sum.eq.0.0 ) s = undef + +! Invert Streamfunction for grads output (in order to be bottom=>top) +! ------------------------------------------------------------------- + do k=1,lm + do j=1,jm + v(j,k) = s(j,lm-k+1) + enddo + enddo + + do k=1,lm + do j=1,jm + s(j,k) = v(j,k) + enddo + enddo + + return + end + subroutine residual ( v0,vpthp0,th0,w0,p0,jm,lm,res,vstar,wstar,wmean,weddy,undef ) + use MAPL_ConstantsMod + implicit none + integer j,k,L,jm,lm + real pi,dp,a,g,H,ps,ts,rhos,z,phi,undef + real airmw,runiv,cpd,rgas,akap, sum + + real v0(jm,lm), v(jm,lm) + real w0(jm,lm), w(jm,lm) + real vpthp0(jm,lm), th0(jm,lm) + real vpthp (jm,lm), th (jm,lm), dthdp(jm,lm) + real stuff (jm,lm) + real res (jm,lm) + real vtlda(jm,lm) + real vstar(jm,lm) + real wstar(jm,lm), wmean(jm,lm), weddy(jm,lm) + real s (jm,lm) + real p0(jm,lm), p(jm,lm), rho0(jm,lm) + real pm(jm,lm) + real cosp(jm), dum(jm,lm) + real ddcosp(jm,lm) + real the(jm,0:lm) + real ple(jm,0:lm) + real stuffe(jm,0:lm) + real delp(jm, lm) + integer method + logical defined + + PARAMETER ( AIRMW = MAPL_AIRMW ) + PARAMETER ( RUNIV = MAPL_RUNIV ) + PARAMETER ( CPD = MAPL_CP ) + PARAMETER ( RGAS = RUNIV/AIRMW ) + PARAMETER ( AKAP = MAPL_KAPPA ) + +! Invert v,th,vpthp, and P level index (in order to be top=>bottom) +! ----------------------------------------------------------------- + do L=1,lm + w(:,L) = w0(:,lm-L+1) + v(:,L) = v0(:,lm-L+1) + p(:,L) = p0(:,lm-L+1) + th(:,L) = th0(:,lm-L+1) + vpthp(:,L) = vpthp0(:,lm-L+1) + enddo + + pi = 4.*atan(1.) + dp = pi/(jm-1) + a = MAPL_RADIUS + g = MAPL_GRAV + H = 7000.0 + ps = 1000.0 + ts = 240.0 + rhos = ps/(rgas*ts) + + Method = 0 +! print *, ' rhos = ', rhos +! print *, '1/rhos = ',1.0/rhos + +! Compute Mean Air Density +! ------------------------ + do L=1,lm + do j=1,jm + z = -H*log(p(j,L)/ps) + rho0(j,L) = rhos*exp(-z/H) + enddo + enddo + + do j=1,jm + phi = -pi/2 + (j-1)*dp + cosp(j) = cos(phi) + enddo + +! Compute Edge Pressures and Thickness +! ------------------------------------ + the(:,0) = th(:,1) + ple(:,0) = max( 0.0, p(:,1) - 0.5*( p(:,2)-p(:,1) ) ) + do L=1,lm-1 + do j=1,jm + ple(j,L) = ( p(j,L)+ p(j,L+1) )*0.5 + enddo + enddo + ple(:,lm) = p(:,lm) + 0.5*( p(:,lm)-p(:,lm-1) ) + + do L=1,lm + delp(:,L) = ple(:,L)-ple(:,L-1) + enddo + +! Compute Mid-Point of PLE Edge Values +! ------------------------------------ + do L=1,lm + pm(:,L) = ( ple(:,L)+ple(:,L-1) )*0.5 + enddo + + +! Compute THE at PLE Edge Values +! ------------------------------ + call map1_cubic( lm,p,th, lm+1,ple,the, jm, Method, undef) + + +! Compute D(Theta)/DZ (with a forced minimum to prevent dthdz => 0) +! ----------------------------------------------------------------- + call compute_d_dp( the,delp,jm,lm,undef,stuff ) + call map1_cubic( lm,pm,stuff, lm,p,dthdp, jm, Method, undef) + + do L=1,lm + do j=1,jm + if( defined(the(j,L ),undef ) .and. & + defined(the(j,L-1),undef ) ) then + dthdp(j,L) = min( -0.003, dthdp(j,L) ) + else + dthdp(j,L) = undef + endif + enddo + enddo + +! Compute Vtlda based on D(rho*vpthp/dthdz)/DZ +! -------------------------------------------- + do L=1,lm + do j=1,jm + if( defined(dthdp(j,L),undef) .and. & + defined(vpthp(j,L),undef) .and. & + dthdp(j,L).ne.0.0 ) then + stuff(j,L) = rho0(j,L)*vpthp(j,L)/(p(j,L)*dthdp(j,L)) + else + stuff(j,L) = undef + endif + enddo + enddo + + call map1_cubic( lm,p,stuff, lm+1,ple,stuffe, jm, Method, undef) + call compute_d_dp( stuffe,delp,jm,lm,undef,stuff ) + call map1_cubic( lm,pm,stuff, lm,p,vtlda, jm, Method, undef) + + do L=1,lm + do j=1,jm + if( defined(vtlda(j,L),undef) ) then + vtlda(j,L) = p(j,L)/rho0(j,L) * vtlda(j,L) + endif + enddo + enddo + +! Compute Vstar +! ------------- + do L=1,lm + do j=1,jm + if( defined( vtlda(j,L),undef) .and. & + defined( v(j,L),undef) ) then + vstar(j,L) = v(j,L) - vtlda(j,L) + else + vstar(j,L) = undef + endif + enddo + enddo + +! Construct Residual Streamfunction from Vstar +! -------------------------------------------- + do k=1,lm + dum(:,1) = 0.0 + do L=1,k + do j=1,jm + if( defined(vstar(j,L),undef) ) then + dum(j,1) = dum(j,1) + vstar(j,L)*delp(j,L)*cosp(j)*rho0(j,L)*H/p(j,L) + endif + enddo + enddo + res(:,k) = dum(:,1) + enddo + + sum = 0.0 + do k=1,lm + do j=1,jm + sum = sum + abs( res(j,k) ) + enddo + enddo + if( sum.eq.0.0 ) res = undef + +! Invert Streamfunction and Vstar for grads output (in order to be bottom=>top) +! ----------------------------------------------------------------------------- + do L=1,lm + do j=1,jm + dum(j,L) = res(j,LM-L+1) + enddo + enddo + do L=1,lm + do j=1,jm + res(j,L) = dum(j,L) + enddo + enddo + + do L=1,lm + do j=1,jm + dum(j,L) = vstar(j,LM-L+1) + enddo + enddo + do L=1,lm + do j=1,jm + vstar(j,L) = dum(j,L) + enddo + enddo + +! Compute D(cos*vpthp/dthdz)/(a*cos*Dphi) +! --------------------------------------- + do L=1,lm + do j=1,jm + if( defined(vpthp(j,L),undef) .and. & + defined(dthdp(j,L),undef) .and. & + dthdp(j,L).ne.0.0 ) then + stuff(j,L) = -H*vpthp(j,L)/(p(j,L)*dthdp(j,L)) + else + stuff(j,L) = undef + endif + enddo + enddo + + call compute_d_dphi( stuff,jm,lm,undef,ddcosp ) + +! Compute Wstar +! ------------- + do L=1,lm + do j=1,jm + if( defined(ddcosp(j,L),undef) ) then + wmean(j,Lm-L+1) = w(j,L) + weddy(j,Lm-L+1) = ddcosp(j,L) + wstar(j,Lm-L+1) = w(j,L) + ddcosp(j,L) + else + wstar(j,Lm-L+1) = undef + wmean(j,Lm-L+1) = undef + weddy(j,Lm-L+1) = undef + endif + enddo + enddo + + + return + end + + subroutine make_w ( v0,p0,jm,lm,w,undef ) + use MAPL_ConstantsMod + implicit none + integer j,k,L,jm,lm + + real v0(jm,lm), v(jm,lm) + real p0(jm,lm), p(jm,lm) + real w(jm,lm), rho0(jm,lm) + real s(jm,lm), cosp(jm) + real dum(jm) + real dvcos_dphi(jm,lm) + real ple(jm,0:lm) + real delp(jm, lm) + logical defined + real airmw,runiv,cpd,rgas,akap + real pi,dp,a,g,H,ps,ts,rhos,phi,z,undef + + PARAMETER ( AIRMW = MAPL_AIRMW ) + PARAMETER ( RUNIV = MAPL_RUNIV ) + PARAMETER ( CPD = MAPL_CP ) + PARAMETER ( RGAS = RUNIV/AIRMW ) + PARAMETER ( AKAP = MAPL_KAPPA ) + +! Invert v and P level index (in order to be top=>bottom) +! ------------------------------------------------------- + do L=1,lm + v(:,L) = v0(:,lm-L+1) + p(:,L) = p0(:,lm-L+1) + enddo + + pi = 4.*atan(1.) + dp = pi/(jm-1) + a = MAPL_RADIUS + g = MAPL_GRAV + H = 7000.0 + ps = 1000.0 + ts = 240.0 + rhos = ps/(rgas*ts) + + do L=1,lm + do j=1,jm + phi = -pi/2 + (j-1)*dp + cosp(j) = cos(phi) + z = -H*log(p(j,L)/ps) + rho0(j,L) = rhos*exp(-z/H) + enddo + enddo + + do L=1,lm + do j=1,jm + if( j.gt.1 .and. j.lt.jm ) then + if( defined(v(j+1,L),undef) .and. & + defined(v(j-1,L),undef) ) then + dvcos_dphi(j,L) = ( v(j+1,L)*cosp(j+1)-v(j-1,L)*cosp(j-1) )/(2*dp) + else + dvcos_dphi(j,L) = undef + endif + else + dvcos_dphi(j,L) = undef + endif + enddo + enddo + +! Compute Edge Pressures and Thickness +! ------------------------------------ + ple(:,0) = max( 0.0, p(:,1) - 0.5*( p(:,2)-p(:,1) ) ) + do L=1,lm-1 + ple(:,L) = ( p(:,L)+ p(:,L+1) )*0.5 + enddo + ple(:,lm) = p(:,lm) + 0.5*( p(:,lm)-p(:,lm-1) ) + + do L=1,lm + delp(:,L) = ple(:,L)-ple(:,L-1) + enddo + +! Construct W from Continuity +! --------------------------- + do k=1,lm + dum(:) = 0.0 + do L=1,k + do j=1,jm + phi = -pi/2+(j-1)*dp + if( dvcos_dphi(j,L).ne.undef ) then + dum(j) = dum(j) + dvcos_dphi(j,L)*delp(j,L)*rho0(j,L)*H/(p(j,L)*a*cosp(j)) + endif + enddo + enddo + s(:,k) = dum(:)/rho0(:,k) + enddo + +! Invert Streamfunction for grads output (in order to be bottom=>top) +! ------------------------------------------------------------------- + do k=1,lm + do j=1,jm + w(j,k) = s(j,lm-k+1) + enddo + enddo + + return + end + + ! ************************************************************************************************************ + + ! function defined ( q,undef ) + ! implicit none + ! logical defined + ! real q,undef + ! defined = abs(q-undef).gt.0.1*undef + ! return + ! end function defined + + ! ************************************************************************************************************ + + subroutine compute_edge( q,p,pe,jm,lm,undef,qe ) + implicit none + integer j,L,jm,lm + real undef + logical defined + real q(jm, lm), p(jm, lm) + real qe(jm,0:lm), pe(jm,0:lm) + + qe(:,0) = q(:,1) + do L=1,lm-1 + do j=1,jm + qe(j,L) = undef + if( defined( q(j,L ),undef) ) qe(j,L) = q(j,L) + if( defined( q(j,L+1),undef) ) qe(j,L) = q(j,L+1) + + ! Linear Interpolation to Pressure Edge + ! ------------------------------------- + if( defined( q(j,L+1),undef) .and. defined( q(j,L),undef) ) then + qe(j,L) = q(j,L) + ( q(j,L+1)-q(j,L) ) * log(pe(j,L)/p(j,L)) / log(p(j,L+1)/p(j,L)) + endif + + enddo + enddo + qe(:,lm) = undef + where( abs( q(:,lm)-undef).gt.0.1 .and. abs( qe(:,lm-1)-undef).gt.0.1 ) + qe(:,lm) = qe(:,lm-1) + ( q(:,lm)- qe(:,lm-1) ) * ( log(pe(:,lm)/pe(:,lm-1)) )/( log(p(:,lm)/pe(:,lm-1)) ) + endwhere + + return + end + + ! ************************************************************************************************************ + + subroutine compute_d_dp( qe,dp,jm,lm,undef,dqdp ) + implicit none + integer j,L,jm,lm + real undef + logical defined + real dp(jm, lm) + real qe(jm,0:lm) + real dqdp(jm,lm) + + do L=1,lm + do j=1,jm + if( defined(qe(j,L-1),undef) .and. defined(qe(j,L),undef) ) then + dqdp(j,L) = ( qe(j,L)-qe(j,L-1) )/ dp(j,L) + else + dqdp(j,L) = undef + endif + enddo + enddo + + return + end + + ! ************************************************************************************************************ + + subroutine compute_d_dphi( q,jm,lm,undef,dqdphi ) + use MAPL_ConstantsMod + implicit none + integer j,L,jm,lm + real q(jm,lm) + real dqdphi(jm,lm) + + real undef, dphi, phi, pi, a + real stuff(jm,lm) + logical defined + + pi = 4.*atan(1.) + dphi = pi/(jm-1) + a = MAPL_RADIUS + + do L=1,lm + do j=1,jm + phi = -pi/2 + (j-1)*dphi + if( defined(q(j,L),undef) ) then + stuff(j,L) = q(j,L)*cos(phi) + else + stuff(j,L) = undef + endif + enddo + enddo + + do L=1,lm + dqdphi(1 ,L) = undef + dqdphi(jm,L) = undef + do j=2,jm-1 + phi = -pi/2 + (j-1)*dphi + if( defined(stuff(j+1,L),undef) .and. & + defined(stuff(j-1,L),undef) ) then + dqdphi(j,L) = ( stuff(j+1,L)-stuff(j-1,L) )/(a*cos(phi)*2*dphi) + else + dqdphi(j,L) = undef + endif + enddo + enddo + + return + end + + ! ************************************************************************************************************ + + subroutine map1_cubic( km, pe1, q1, kn, pe2, q2, jm, Method, undef) + use MAPL_ConstantsMod + implicit none + + real, intent(in) :: undef + integer, intent(in) :: Method ! 0: Linear in P + ! 1: Linear in Log(P) + ! 2: Linear in P**kappa + integer, intent(in) :: jm ! Latitude dimension + integer, intent(in) :: km ! Original vertical dimension + integer, intent(in) :: kn ! Target vertical dimension + + real, intent(in) :: pe1(jm,km) ! pressure at mid-layers + ! in the original vertical coordinate + + real, intent(in) :: pe2(jm,kn) ! pressure at mid-layers + ! in the new vertical coordinate + + real, intent(in) :: q1(jm,km) ! Field input + real, intent(inout):: q2(jm,kn) ! Field output + +! !DESCRIPTION: +! +! Perform Cubic Interpolation in the vertical +! ------------------------------------------- +! pe1: pressure associated with q1 +! pe2: pressure associated with q2 +! +!----------------------------------------------------------------------- +! + real qx(jm,km) + real logpl1(jm,km) + real logpl2(jm,kn) + real dlogp1(jm,km) + real am2,am1,ap0,ap1,P,PLP1,PLP0,PLM1,PLM2,DLP0,DLM1,DLM2 + + integer j, k, LM2,LM1,LP0,LP1 + logical defined + + real airmw,runiv,cpd,rgas,akap + PARAMETER ( AIRMW = MAPL_AIRMW ) + PARAMETER ( RUNIV = MAPL_RUNIV ) + PARAMETER ( CPD = MAPL_CP ) + PARAMETER ( RGAS = RUNIV/AIRMW ) + PARAMETER ( AKAP = MAPL_KAPPA ) + +! Initialization +! -------------- + + select case (Method) + + ! Linear in P + ! ----------- + case(0) + do k=1,km + qx(:,k) = q1(:,k) + logpl1(:,k) = pe1(:,k) + enddo + do k=1,kn + logpl2(:,k) = pe2(:,k) + enddo + + do k=1,km-1 + dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k) + enddo + + ! Linear in Log(P) + ! ---------------- + case(1) + do k=1,km + qx(:,k) = q1(:,k) + logpl1(:,k) = log( pe1(:,k) ) + enddo + do k=1,kn + logpl2(:,k) = log( pe2(:,k) ) + enddo + + do k=1,km-1 + dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k) + enddo + + ! Linear in P**kappa + ! ------------------ + case(2) + do k=1,km + qx(:,k) = q1(:,k) + logpl1(:,k) = exp( akap*log( pe1(:,k) )) + enddo + do k=1,kn + logpl2(:,k) = exp( akap*log( pe2(:,k) )) + enddo + + do k=1,km-1 + dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k) + enddo + + end select + +! Interpolate Q1 onto target Pressures +! ------------------------------------ + do j=1,jm + do k=1,kn + LM1 = 1 + LP0 = 1 + do while( LP0.le.km ) + if (logpl1(j,LP0).lt.logpl2(j,k)) then + LP0 = LP0+1 + else + exit + endif + enddo + LM1 = max(LP0-1,1) + LP0 = min(LP0, km) + +! Extrapolate Linearly above first model level +! -------------------------------------------- + if( LM1.eq.1 .and. LP0.eq.1 ) then + q2(j,k) = qx(j,1) + if( defined(qx(j,2),undef) ) q2(j,k) = qx(j,2) + if( defined(qx(j,1),undef) .and. defined(qx(j,2),undef) ) then + q2(j,k) = qx(j,1) + ( qx(j,2)-qx(j,1) )*( logpl2(j,k)-logpl1(j,1) ) & + /( logpl1(j,2)-logpl1(j,1) ) + endif + +! Extrapolate Linearly below last model level +! ------------------------------------------- + else if( LM1.eq.km .and. LP0.eq.km ) then + q2(j,k) = qx(j,km) + if( defined(qx(j,km-1),undef) ) q2(j,k) = qx(j,km-1) + if( defined(qx(j,km-1),undef) .and. defined(qx(j,km),undef) ) then + q2(j,k) = qx(j,km) + ( qx(j,km)-qx(j,km-1) )*( logpl2(j,k )-logpl1(j,km ) ) & + /( logpl1(j,km)-logpl1(j,km-1) ) + endif + +! Interpolate Linearly between levels 1 => 2 and km-1 => km +! --------------------------------------------------------- + else if( LM1.eq.1 .or. LP0.eq.km ) then + q2(j,k) = qx(j,LP0) + if( defined(qx(j,LM1),undef) ) q2(j,k) = qx(j,LM1) + if( defined(qx(j,LP0),undef) .and. defined(qx(j,LM1),undef) ) then + q2(j,k) = qx(j,LP0) + ( qx(j,LM1)-qx(j,LP0) )*( logpl2(j,k )-logpl1(j,LP0) ) & + /( logpl1(j,LM1)-logpl1(j,LP0) ) + endif + +! Interpolate Cubicly between other model levels +! ---------------------------------------------- + else + LP1 = LP0+1 + LM2 = LM1-1 + P = logpl2(j,k) + PLP1 = logpl1(j,LP1) + PLP0 = logpl1(j,LP0) + PLM1 = logpl1(j,LM1) + PLM2 = logpl1(j,LM2) + DLP0 = dlogp1(j,LP0) + DLM1 = dlogp1(j,LM1) + DLM2 = dlogp1(j,LM2) + + ap1 = (P-PLP0)*(P-PLM1)*(P-PLM2)/( DLP0*(DLP0+DLM1)*(DLP0+DLM1+DLM2) ) + ap0 = (PLP1-P)*(P-PLM1)*(P-PLM2)/( DLP0* DLM1 *( DLM1+DLM2) ) + am1 = (PLP1-P)*(PLP0-P)*(P-PLM2)/( DLM1* DLM2 *(DLP0+DLM1 ) ) + am2 = (PLP1-P)*(PLP0-P)*(PLM1-P)/( DLM2*(DLM1+DLM2)*(DLP0+DLM1+DLM2) ) + + if( defined(qx(j,LP1),undef) .and. & + defined(qx(j,LP0),undef) .and. & + defined(qx(j,LM1),undef) .and. & + defined(qx(j,LM2),undef) ) then + q2(j,k) = ap1*qx(j,LP1) + ap0*qx(j,LP0) + am1*qx(j,LM1) + am2*qx(j,LM2) + + else if( defined(qx(j,LP0),undef) .and. defined(qx(j,LM1),undef) ) then + q2(j,k) = qx(j,LP0) + ( qx(j,LM1)-qx(j,LP0) )*( logpl2(j,k )-logpl1(j,LP0) ) & + /( logpl1(j,LM1)-logpl1(j,LP0) ) + + else + q2(j,k) = undef + endif + + endif + + enddo + enddo + + return + end subroutine map1_cubic diff --git a/post/gcmclim.script b/post/gcmclim.script index c10b2d0b..be6be94e 100755 --- a/post/gcmclim.script +++ b/post/gcmclim.script @@ -106,7 +106,7 @@ endif echo Use: NCPUS = $NCPUS if( $NCPUS == "NULL" | $NCPUS == 1 ) then - set TIMEAVE_x = "$GEOSUTIL/bin/time_ave.x" + set TIMEAVE_x = "$GEOSUTIL/bin/time_ave.x -NO_TEM true" else if ( ! $?RUN_CMD ) then echo " " @@ -114,7 +114,7 @@ echo Use: NCPUS = $NCPUS echo " " exit 1 endif - set TIMEAVE_x = "$RUN_CMD $NCPUS $GEOSUTIL/bin/time_ave.x" + set TIMEAVE_x = "$RUN_CMD $NCPUS $GEOSUTIL/bin/time_ave.x -NO_TEM true" endif diff --git a/post/gcmpost.script b/post/gcmpost.script index 64087fa1..60db4576 100755 --- a/post/gcmpost.script +++ b/post/gcmpost.script @@ -679,13 +679,15 @@ sleep 5 echo " Performing time average ..." if(! -e time_ave.rc ) ln -s $GEOSUTIL/post/time_ave.rc . if( $NCPUS == "NULL" | $NCPUS == 1 | !($?RUN_CMD) ) then - $GEOSUTIL/bin/time_ave_nompi.x -hdf $locals -rc time_ave.rc \ - -tag $expid.$collection.monthly \ - -d $expid.$collection.diurnal $STRICT $IGNORE_NAN + $GEOSUTIL/bin/time_ave_nompi.x -hdf $locals -rc time_ave.rc \ + -tag $expid.$collection.monthly \ + -d $expid.$collection.diurnal \ + -expid $expid $STRICT $IGNORE_NAN else - $RUN_CMD $NCPUS $GEOSUTIL/bin/time_ave.x -hdf $locals -rc time_ave.rc \ - -tag $expid.$collection.monthly \ - -d $expid.$collection.diurnal $STRICT $IGNORE_NAN + $RUN_CMD $NCPUS $GEOSUTIL/bin/time_ave.x -hdf $locals -rc time_ave.rc \ + -tag $expid.$collection.monthly \ + -d $expid.$collection.diurnal \ + -expid $expid $STRICT $IGNORE_NAN endif endif @@ -743,19 +745,20 @@ foreach collection ( $collections ) if( $streamdates[1] != 0 ) then foreach date ( $streamdates ) - if( -e $SOURCE/holding/$collection/$date ) then + if( -e $SOURCE/holding/$collection/$date ) then set locdir = $SOURCE/holding/$collection/$date cd $SOURCE/holding/$collection/$date - set archfile = $SOURCE/archive/archive_commands.$collection.$date.`date +%Y%m%d_%H%M%S` + set archfile = $SOURCE/archive/archive_commands.$collection.$date.`date +%Y%m%d_%H%M%S` /bin/rm -f $archfile - touch $archfile + touch $archfile if( ${#monthly_attribute} == 0 ) then - set num = `/bin/ls -1 $expid.$collection.monthly.* | wc -l` + set num = `/bin/ls -1 $expid.$collection.monthly.* | wc -l` else - set num = `/bin/ls -1 $expid.$collection.* | wc -l` + set num = `/bin/ls -1 $expid.$collection.* | wc -l` endif + set TEMnum = `/bin/ls -1 $expid.TEM_Diag.monthly.* | wc -l` if( $num != 0 ) then @@ -951,6 +954,75 @@ foreach collection ( $collections ) end endif + + if( $TEMnum != 0 ) then + + mkdir -p $SOURCE/TEM_Diag + + if( $check_plot == 'NULL' ) set check_plot = TRUE + echo "set success = TRUE" >> $archfile + echo "if( -e $SOURCE/holding/$collection/$date ) then " >> $archfile + echo " cd $SOURCE/holding/$collection/$date" >> $archfile + echo "wait" >> $archfile + + # Archive TEM_Diag Monthlies + # -------------------------- + @ monthly_date_node = $date_node + 1 # Since monthly is simply appended to collection name + + set locals = `/bin/ls -1 | grep "${expid}.TEM_Diag.*${ext}" | grep monthly | grep -v diurnal | grep -v partial` + + echo "Archiving Monthlies for Collection: TEM_Diag for Date: $date" + + foreach local ($locals) + set fdate = `echo $local | cut -d "." -f${monthly_date_node}` # Assumed date format: YYYYMMDD_HHMMSS + set year = `echo $fdate | cut -c1-4` + set month = `echo $fdate | cut -c5-6` + set day = `echo $fdate | cut -c7-8` + set hour = `echo $fdate | cut -c10-11` + + set archive = `echo $archives[$m] | sed -e "s/%c/TEM_Diag/ g" | \ + sed -e "s/%y4/$year/ g" | \ + sed -e "s/%m2/$month/ g" | \ + sed -e "s/%d2/$day/ g" | \ + sed -e "s/%h2/$hour/ g"` + + $cpvar $local $SOURCE/TEM_Diag + + echo "if( -e $local ) then" >> $archfile + echo "ssh ${MASTOR} mkdir -p ${MASDIR}/$archive" >> $archfile + echo 'echo " "' >> $archfile + echo "echo 'Archiving: '$local" >> $archfile + + if( $HOST == 'pleiades' ) then + echo "ssh ${MASTOR} $scpvar $locdir/$local ${MASDIR}/$archive" >> $archfile + + echo 'if ($status == 0) then ' >> $archfile + echo 'echo " Archive ... PASS "' >> $archfile + if( $FSEGMENT == 00000000 ) echo "/bin/rm -f $local" >> $archfile + echo 'else' >> $archfile + echo 'echo " Archive ... FAIL!"' >> $archfile + echo 'set success = FALSE' >> $archfile + echo 'endif' >> $archfile + + else if( $HOST == 'discover' ) then + echo "$scpvar $local ${MASTOR}:${MASDIR}/$archive | set archive_status = "'`sed -e "s/ / /g"`' >> $archfile + + echo 'if ( .$archive_status[4] == .OK ) then' >> $archfile + echo 'echo " Archive ... PASS "' >> $archfile + if( $FSEGMENT == 00000000 ) echo "/bin/rm -f $local" >> $archfile + echo 'else' >> $archfile + echo 'echo " Archive ... FAIL!"' >> $archfile + echo 'set success = FALSE' >> $archfile + echo 'endif' >> $archfile + + echo 'unset archive_status' >> $archfile + endif + echo 'endif' >> $archfile + + end + + endif + # Archive Restarts (on First Stream Only) # --------------------------------------- if( $collection == $streams[1] ) then @@ -1141,6 +1213,90 @@ foreach collection ( $collections ) endif end # End ForEach Collection Loop + # Create DDF File for TEM_Diags + # ----------------------------- + if(-e $SOURCE/TEM_Diag ) then + cd $SOURCE/TEM_Diag + set TEM_Diags = `/bin/ls -1 $expid.TEM_Diag.monthly* | grep -v clim | grep -v partial | grep -v ddf` + set num = `/bin/ls -1 $TEM_Diags | wc -l` + + echo " " + echo "TEM_Diag files: $TEM_Diags" + echo "Number: $num" + echo " " + + if( $num != 0 ) then + + # Extract Extension (last node) from Monthly Mean File + # ---------------------------------------------------- + @ m = 1 + while( $m != 0 ) + set dum = `echo $TEM_Diags[1] | cut -d "." -f$m` + if( .$dum != . ) then + set ext = $dum + @ m = $m + 1 + else + @ m = 0 + endif + end + echo "Extension: $ext" + echo " " + + set edate = `echo $TEM_Diags[$num] | cut -d "." -f${date_node}` + set eyear = `echo $edate | cut -c1-4` + set emonth = `echo $edate | cut -c5-6` + set eday = `echo $edate | cut -c7-8` + set ehour = `echo $edate | cut -c10-11` + + echo " edate: $edate" + echo " eyear: $eyear" + echo "emonth: $emonth" + echo " eday: $eday" + echo " ehour: $ehour" + + set date = `echo $TEM_Diags[1] | cut -d "." -f${date_node}` + set year = `echo $date | cut -c1-4` + set month = `echo $date | cut -c5-6` + set day = `echo $date | cut -c7-8` + set hour = `echo $date | cut -c10-11` + + echo " date: $date" + echo " year: $year" + echo " month: $month" + echo " day: $day" + echo " hour: $hour" + + set expdsc = `grep EXPDSC: $HISTORYRC | cut -d" " -f2` + set timinc = "mo" + + if( $month == "01" ) set MON = JAN + if( $month == "02" ) set MON = FEB + if( $month == "03" ) set MON = MAR + if( $month == "04" ) set MON = APR + if( $month == "05" ) set MON = MAY + if( $month == "06" ) set MON = JUN + if( $month == "07" ) set MON = JUL + if( $month == "08" ) set MON = AUG + if( $month == "09" ) set MON = SEP + if( $month == "10" ) set MON = OCT + if( $month == "11" ) set MON = NOV + if( $month == "12" ) set MON = DEC + + @ nmonths = ( $eyear - $year ) * 12 + ( $emonth - $month ) + 1 + + /bin/rm -f xdf.tabl + echo DSET $SOURCE/TEM_Diag/$expid.TEM_Diag.monthly.%y4%m2.$ext > xdf.tabl + echo TITLE $expdsc >> xdf.tabl + echo OPTIONS template >> xdf.tabl + echo TDEF time $nmonths LINEAR 00:00Z01$MON$year 1$timinc >> xdf.tabl + /bin/ln -sf xdf.tabl ${expid}.TEM_Diag.monthly.ddf + + echo "xdf.tabl: " + cat xdf.tabl + + endif # Endif for TEM_Diags NUM = 0 check + endif # Endif for TEM_Diag Directory check + endif # End PLOTS Option ####################################################################### diff --git a/post/time_ave.F b/post/time_ave.F index 5fb90ffa..d65ca83f 100644 --- a/post/time_ave.F +++ b/post/time_ave.F @@ -107,6 +107,39 @@ program time_ave integer nalias logical, allocatable :: lzstar(:) +C ********************************************************************** +C **** Variables for TEM Diagnostics **** +C ********************************************************************** + + integer, parameter :: TEM_NVARS = 7 + character*256 TEM_NAMES(TEM_NVARS) + logical TEM_FOUND(TEM_NVARS) + integer TEM_BEG (TEM_NVARS) + integer TEM_END (TEM_NVARS) + + real, allocatable :: TEM_VARS(:,:,:,:) + + data TEM_NAMES / 'U','V','T','OMEGA','USVS','USWS','VSTS' / + data TEM_FOUND / '.F.' , + . '.F.' , + . '.F.' , + . '.F.' , + . '.F.' , + . '.F.' , + . '.F.' / + + logical TEM_WRITE + logical NO_TEM + + integer mem_unit,loc,ios + logical lopen + character*256 cbuff1 + character*1 cbuff2(256) + equivalence ( cbuff1,cbuff2 ) + character*256 expid + +C ********************************************************************** + integer NSECF, ntmin, ntcrit, nhmsf, nc NSECF(N) = N/10000*3600 + MOD(N,10000)/100* 60 + MOD(N,100) @@ -156,6 +189,8 @@ program time_ave diurnal = .FALSE. mdiurnal = .FALSE. ignore_nan = .FALSE. + NO_TEM = .FALSE. + expid = 'NULL' allocate ( arg(nargs) ) do n=1,nargs call getarg(n,arg(n)) @@ -171,8 +206,10 @@ program time_ave if( trim(arg(n)).eq.'-ntmin' ) read ( arg(n+1),* ) ntmin if( trim(arg(n)).eq.'-ntod' ) read ( arg(n+1),* ) ntod if( trim(arg(n)).eq.'-ndt' ) read ( arg(n+1),* ) ndt + if( trim(arg(n)).eq.'-expid' ) read ( arg(n+1),* ) expid if( trim(arg(n)).eq.'-strict' ) read ( arg(n+1),* ) strict if( trim(arg(n)).eq.'-noquad' ) lquad = .FALSE. + if( trim(arg(n)).eq.'-NO_TEM' ) NO_TEM = .TRUE. if( trim(arg(n)).eq.'-ignore_nan' ) ignore_nan = .TRUE. if( trim(arg(n)).eq.'-dv'.or. trim(arg(n)).eq.'-mdv') ldquad = .true. @@ -533,7 +570,8 @@ program time_ave print * print *, 'Files: ' do n=1,nfiles - print *, n,trim(fname(n)) + write(6,7004) n,trim(fname(n)) + 7004 format(1x,i4,3x,a) enddo print * if( ntod.eq.-999 ) then @@ -562,10 +600,10 @@ program time_ave k = 0 do n=1,nfiles - if( root ) then - call gfio_open ( trim(fname(n)),1,ID,rc ) - call gfio_diminquireCF ( id,imglobal,jmglobal,lm,ntime,nvars,ngatts,twoDimLat,rc ) - endif + if( root ) then + call gfio_open ( trim(fname(n)),1,ID,rc ) + call gfio_diminquireCF ( id,imglobal,jmglobal,lm,ntime,nvars,ngatts,twoDimLat,rc ) + endif #ifdef mpi call mpi_bcast ( imglobal,1,mpi_integer,0,comm,ierror ) @@ -861,6 +899,17 @@ program time_ave 3001 format(1x,'Writing ',a,' at location ',i6,' into File: ',a) dum(:,:,1:kend) = q(:,:,nloc(n):nloc(n)+kend-1,0) call mpi_gfio_putVar ( id,trim(vname2(n)),nymd0,nhms0,im,jm,kbeg,kend,dum,lattice ) + +c Tests for TEM Diagnostic Variables +c ---------------------------------- + do k=1,TEM_NVARS + if( trim(vname2(n)).eq.TEM_NAMES(k) ) then + TEM_FOUND(k) = .TRUE. + TEM_BEG (k) = nloc(n) + TEM_END (k) = nloc(n)+kend-1 + endif + enddo + enddo c Quadratics @@ -888,8 +937,19 @@ program time_ave endwhere else dum(:,:,1:kend) = q(:,:,nloc(mv):nloc(mv)+kend-1,0) - endif + endif call mpi_gfio_putVar ( id,trim(vname2(mv)),nymd0,nhms0,im,jm,kbeg,kend,dum,lattice ) + +c Tests for TEM Diagnostic Variables +c ---------------------------------- + do k=1,TEM_NVARS + if( trim(vname2(mv)).eq.TEM_NAMES(k) ) then + TEM_FOUND(k) = .TRUE. + TEM_BEG (k) = nloc(mv) + TEM_END (k) = nloc(mv)+kend-1 + endif + enddo + endif enddo @@ -897,10 +957,66 @@ program time_ave call gfio_close ( id,rc ) print * print *, 'Created: ',trim(hdfile) + print *, '-------- ' print * endif call timeend(' Write_AVE') +C ********************************************************************** +C **** Write TEM Diagnostics File **** +C ********************************************************************** + + TEM_WRITE = .not.NO_TEM + do k=1,TEM_NVARS + TEM_WRITE = TEM_WRITE .and. TEM_FOUND(k) + enddo + if( TEM_WRITE ) then + + allocate( TEM_VARS(imglobal,jmglobal,lm,TEM_NVARS) ) + do k=1,TEM_NVARS + do L=TEM_BEG(k),TEM_END(k) + dum(:,:,L-TEM_BEG(k)+1) = q(:,:,L,0) + call gather_2d( TEM_VARS(1,1,L-TEM_BEG(k)+1,k),dum(1,1,L-TEM_BEG(k)+1),lattice ) + enddo + enddo + + if( root ) then + + if( trim(expid).eq.'NULL' ) then + cbuff1 = '/usr/bin/basename ' // trim(fname(1)) // ' > basename.output' + call execute_command_line ( trim(cbuff1) ) + + i = 1 + mem_unit = 0 + dowhile( mem_unit == 0 ) + if( i /= 5 .and. i /= 6 .and. i <= 99 ) then + inquire( unit = i, opened = lopen, iostat = ios ) + if( ios == 0 ) then + if( .not. lopen ) mem_unit = i + endif + endif + i = i+1 + enddo + + open (mem_unit,file='basename.output',form='formatted',access='sequential') + read (mem_unit,'(a)') cbuff1 + close(mem_unit) + + expid = '' + loc = 1 + char = cbuff2(loc) + dowhile (char.ne.'.') + expid = trim(expid) // char + loc = loc+1 + char = cbuff2(loc) + enddo + endif + + call TEM( TEM_VARS,imglobal,jmglobal,lm,TEM_NVARS, + . EXPID,lon,lat,lev,levunits,nymd0,nhms0,timinc,undef ) + endif + endif + C ********************************************************************** C **** Write HDF Monthly Diurnal Output File **** C ********************************************************************** @@ -985,6 +1101,8 @@ program time_ave if( root .and. mdiurnal ) then call gfio_close ( id,rc ) print *, 'Created: ',trim(hdfile) + print *, '-------- ' + print * endif call tick (nymd0,nhms0,ndt) enddo @@ -992,6 +1110,8 @@ program time_ave if( root .and. diurnal ) then call gfio_close ( id,rc ) print *, 'Created: ',trim(hdfile) + print *, '-------- ' + print * endif if( root ) print * From 7500ceec790f9f83396cc7b7b4eb4164d1640ce9 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 10 Jun 2024 15:14:03 -0400 Subject: [PATCH 2/3] Update changelog for 2.1.0 release --- CHANGELOG.md | 17 +++++++++++------ post/gcmpost_CPLFCSTfull.script | 4 ++-- post/gcmpost_CPLFCSTpart.script | 4 ++-- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e0af4a4a..9e6abdcc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,20 +9,25 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Added new BCS version v12 - ### Changed -- Move to use `time_ave_util.x` from MAPL instead of `time_ave.x` - ### Fixed ### Removed -- mask from `read_reynolds.F90` - ### Deprecated +## [2.1.0] - 2024-06-10 + +### Added + +- Added new BCS version v12 +- Updates for TEM diagnostics + +### Removed + +- Removed mask from `read_reynolds.F90` + ## [2.0.8] - 2024-03-29 ### Fixed diff --git a/post/gcmpost_CPLFCSTfull.script b/post/gcmpost_CPLFCSTfull.script index 9246fa9c..f75d9951 100755 --- a/post/gcmpost_CPLFCSTfull.script +++ b/post/gcmpost_CPLFCSTfull.script @@ -583,11 +583,11 @@ foreach collection ( $collections ) echo " Performing time average ..." if(! -e time_ave.rc ) ln -s $GEOSUTIL/post/time_ave.rc . if( $NCPUS == "NULL" | $NCPUS == 1 | !($?RUN_CMD) ) then - $GEOSUTIL/post/time_ave_util.x -hdf $locals -rc time_ave.rc \ + $GEOSUTIL/post/time_ave.x -hdf $locals -rc time_ave.rc \ -tag $expid.$collection.monthly \ -d $expid.$collection.diurnal $STRICT else - $RUN_CMD $NCPUS $GEOSUTIL/post/time_ave_util.x -hdf $locals -rc time_ave.rc \ + $RUN_CMD $NCPUS $GEOSUTIL/post/time_ave.x -hdf $locals -rc time_ave.rc \ -tag $expid.$collection.monthly \ -d $expid.$collection.diurnal $STRICT endif diff --git a/post/gcmpost_CPLFCSTpart.script b/post/gcmpost_CPLFCSTpart.script index 1103e7c2..d5a29027 100755 --- a/post/gcmpost_CPLFCSTpart.script +++ b/post/gcmpost_CPLFCSTpart.script @@ -583,12 +583,12 @@ foreach collection ( $collections ) echo " Performing partial time average ..." if(! -e time_ave.rc ) ln -s $GEOSUTIL/post/time_ave.rc . if( $NCPUS == "NULL" | $NCPUS == 1 | !($?RUN_CMD) ) then - $GEOSUTIL/post/time_ave_util.x -hdf $locals -rc time_ave.rc \ + $GEOSUTIL/post/time_ave.x -hdf $locals -rc time_ave.rc \ -tag $expid.$collection.monthly \ -ntmin 1 \ -d $expid.$collection.diurnal $STRICT else - $RUN_CMD $NCPUS $GEOSUTIL/post/time_ave_util.x -hdf $locals -rc time_ave.rc \ + $RUN_CMD $NCPUS $GEOSUTIL/post/time_ave.x -hdf $locals -rc time_ave.rc \ -tag $expid.$collection.monthly \ -ntmin 1 \ -d $expid.$collection.diurnal $STRICT From 5854160d61895dc7ebd11a108eba067f204972e3 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 10 Jun 2024 15:40:13 -0400 Subject: [PATCH 3/3] Remove hardcoded paths --- plots/grads_util/Create_TEM_ddf.csh | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/plots/grads_util/Create_TEM_ddf.csh b/plots/grads_util/Create_TEM_ddf.csh index c588923a..01d0da03 100755 --- a/plots/grads_util/Create_TEM_ddf.csh +++ b/plots/grads_util/Create_TEM_ddf.csh @@ -2,10 +2,8 @@ setenv ARCH `uname` -setenv SITE NCCS -setenv GEOSDIR /discover/nobackup/ltakacs/TAGS/Jason-4_0_p2/GEOSagcm -setenv GEOSBIN /discover/nobackup/ltakacs/TAGS/Jason-4_0_p2/GEOSagcm/Linux/bin -setenv GEOSUTIL /discover/nobackup/ltakacs/TAGS/Jason-4_0_p2/GEOSagcm/src/GMAO_Shared/GEOS_Util +# NOTE The script below requires GEOSBIN. Usually this should be set in scripts calling +# this but if not, you should set it to the install/bin directory of your model source $GEOSBIN/g5_modules setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:${BASEDIR}/${ARCH}/lib