<![CDATA[ Subroutine PUTCDF(TMIDDLE)
*
* Store selected POM variables in NetCDF file. This routine should
* be called whenever you would normally output to the "print" file.
*
* -------------------------------------------------------------------
* Here's how to build NetCDF output into the POM code:
*
* (1) Add the following line to comblk97.h:
* COMMON/GRIDSTUFF/ XGRID(IM,JM),YGRID(IM,JM), ANG(IM,JM)
*
* (2) Define the centers and angles of the grid cells in pom.f:
* XGRID(I,J) x coordinates
* of the grid cell centers in either projected
* or geographic coordinates (x in km, or
* longitude in decimal degrees).
* YGRID(I,J) y coordinates
* of the grid cell centers in either projected
* or geographic coordinates (y in km, or
* latitude in decimal degrees).
* ANG(I,J) angle in degrees between east and the "i" grid
* direction, measured in a counter-clockwise direction
*
* (3) Add this routine (PUTCDF.F) into the POM code.
*
* (4) Call PUTCDF just after the line in pom97.f that reads
* just before the PRXY calls in loop 9000
*
* call putcdf(time)
*
* (5) Get the NetCDF library from
* ftp://unidata.ucar.edu/pub/netcdf
* and build it!
*
* (6) Copy the netcdf.inc library to the POM source code
*
* (7) Link to the netcdf library when you compile POM
* e.g. f77 -O -o pom97 pom97.f /usr/local/lib/libnetcdf.a
*
*---------------------------------------------------------------
* Rich Signell | rsignell@crusty.er.usgs.gov
* U.S. Geological Survey | (508) 457-2229 | FAX (508) 457-2310
* 384 Woods Hole Road | http://crusty.er.usgs.gov/rsignell.html
* Woods Hole, MA 02543-1598 |
*
* Last Modified: Fri Sep 24 14:20:18 EDT 1999
Include 'netcdf.inc'
Include 'comblk97.h'
Real KMMIN, KMMAX, KMSC, KMOF
Character*80 COM1
*
* Many variables have a relatively small dynamic range, and are
* therefore scaled and stored as short integers to save space.
* The scaling is determined by the specified MIN and MAX of
* these variables, which is defined in the following lines.
* These may need to be changed for your particular case!
*
* elevation range: -4 to 4 meters
Parameter (ELMIN = -4.,ELMAX = 4.)
* salinity range -5 to 40 PSU
Parameter (SALMIN = -5,SALMAX = 40.)
* temperature range -10 to 40 C
Parameter (TMPMIN = -10.,TMPMAX = 40.)
* u and v components of velocity -3 to 3 m/s
Parameter (VELMIN = -3.,VELMAX = 3.)
* w component of velocity -.002 m/s to .002 m/s
Parameter (WMIN = -.002,WMAX = 0.002)
Parameter (TAUMIN=-10,TAUMAX=10)
Parameter (KMMIN = 0.,KMMAX = 0.2)
Parameter (AMMIN = 0.,AMMAX = 1000.)
Parameter (HMIN = 0.,HMAX = 280.)
Parameter (CMIN = 0., CMAX = 40.)
Parameter (HFMIN = -4000.,HFMAX = 4000.)
*
Integer IRET, BASE_DATE(3)
* netCDF id
Integer CDFID
* dimension ids
Integer XPOSDIM, YPOSDIM, ZPOSDIM, TIMEDIM
* variable ids
Integer XID, YID, SIGMAID, DEPTHID, TID, ELID, UID, VID, WID,
* SALID, TMPID, ANGID, H1ID, H2ID, TXID,
* TYID, KMID, AMID, HFID, CCID,CDID
* variable shapes
Integer DIMS(4)
* corners and edge lengths
Integer CORNER(4), EDGES(4)
* attribute vectors
Real FLOATVAL(1), VRNG(2), htmp(im,jm)
* short integer workspace arrays
Integer*2 WS0(IM,JM)
Integer*2 WS1(IM,JM,KB)
* count the number of times this routine is called
Data ICDF /0/
Character*1 VSX, VSY, PTU, PTV, PTW, PTAM, PTS, PTT, PRHO, PTQ2,
* PTL, PTKM, PTUF, PTVF, PTHF, PTC, PTCD
Save
If (ICDF.EQ.0) Then
* Example of simple descriptor for run -- something like this would
* normally be supplied by an input file
COM1='Demo POM run'
* If you see a variable listed below, you can choose not to save it
* in the output file by changing the 'y' to 'n'
PTU = 'y'
PTV = 'y'
PTW = 'y'
PTS = 'y'
PTT = 'y'
PTW = 'y'
PTAM = 'y'
PTKM = 'y'
* calculate scale factors and offsets based on expected range of data
ELOF = (ELMAX+ELMIN) * .5
HFOF = (HFMAX+HFMIN) * .5
SALOF = (SALMAX+SALMIN) * .5
TMPOF = (TMPMAX+TMPMIN) * .5
COF = (CMAX+CMIN) * .5
VELOF = (VELMAX+VELMIN) * .5
WOF = (WMAX+WMIN) * .5
KMOF = (KMMAX+KMMIN) * .5
AMOF = (AMMAX+AMMIN) * .5
RANGE = 32767 - (-32767)
ELSC = (ELMAX-ELMIN) / RANGE
HFSC = (HFMAX-HFMIN) / RANGE
SALSC = (SALMAX-SALMIN) / RANGE
TMPSC = (TMPMAX-TMPMIN) / RANGE
CSC = (CMAX-CMIN) / RANGE
VELSC = (VELMAX-VELMIN) / RANGE
WSC = (WMAX-WMIN) / RANGE
KMSC = (KMMAX-KMMIN) / RANGE
AMSC = (AMMAX-AMMIN) / RANGE
TAUSC = (TAUMAX-TAUMIN) / RANGE
TAUOF = (TAUMAX+TAUMIN) * 0.5
* hard-wire: save heat flux
PTHF = 'y'
* hard-wire save concentration
PTC = 'y'
* enter define mode
CDFID = NCCRE('pom97.cdf',NCCLOB,IRET)
IFILL = NCSFIL(CDFID,NCNOFILL,IRET)
* define dimensions
XPOSDIM = NCDDEF(CDFID,'xpos',IM,IRET)
YPOSDIM = NCDDEF(CDFID,'ypos',JM,IRET)
ZPOSDIM = NCDDEF(CDFID,'zpos',KB,IRET)
TIMEDIM = NCDDEF(CDFID,'time',NCUNLIM,IRET)
* define variables
* 1D Vars
DIMS(1) = TIMEDIM
TID = NCVDEF(CDFID,'time',NCFLOAT,1,DIMS,IRET)
DIMS(1) = ZPOSDIM
SIGMAID = NCVDEF(CDFID,'sigma',NCFLOAT,1,DIMS,IRET)
* 2D Vars
DIMS(2) = YPOSDIM
DIMS(1) = XPOSDIM
XID = NCVDEF(CDFID,'x',NCFLOAT,2,DIMS,IRET)
YID = NCVDEF(CDFID,'y',NCFLOAT,2,DIMS,IRET)
H1ID = NCVDEF(CDFID,'h1',NCFLOAT,2,DIMS,IRET)
H2ID = NCVDEF(CDFID,'h2',NCFLOAT,2,DIMS,IRET)
DEPTHID = NCVDEF(CDFID,'depth',NCFLOAT,2,DIMS,IRET)
ANGID = NCVDEF(CDFID,'ang',NCFLOAT,2,DIMS,IRET)
* 3D Vars
DIMS(3) = TIMEDIM
DIMS(2) = YPOSDIM
DIMS(1) = XPOSDIM
ELID = NCVDEF(CDFID,'elev',NCSHORT,3,DIMS,IRET)
* 1D, 2D, and 3D attributes
Call NCAPTC(CDFID,XID,'long_name',NCCHAR,7,'easting',IRET)
Call NCAPTC(CDFID,XID,'units',NCCHAR,6,'meters',IRET)
Call NCAPTC(CDFID,H1ID,'long_name',NCCHAR,8,'x metric',IRET)
Call NCAPTC(CDFID,H1ID,'units',NCCHAR,6,'meters',IRET)
Call NCAPTC(CDFID,YID,'long_name',NCCHAR,8,'northing',IRET)
Call NCAPTC(CDFID,YID,'units',NCCHAR,6,'meters',IRET)
Call NCAPTC(CDFID,H2ID,'long_name',NCCHAR,8,'y metric',IRET)
Call NCAPTC(CDFID,H2ID,'units',NCCHAR,6,'meters',IRET)
Call NCAPTC(CDFID,SIGMAID,'long_name',NCCHAR,36,
* 'stretched vertical coordinate levels',IRET)
Call NCAPTC(CDFID,SIGMAID,'units',NCCHAR,17,'fraction of depth',
* IRET)
Call NCAPTC(CDFID,DEPTHID,'long_name',NCCHAR,10,'Bathymetry',
* IRET)
Call NCAPTC(CDFID,DEPTHID,'units',NCCHAR,6,'meters',IRET)
Call NCAPTC(CDFID,ANGID,'long_name',NCCHAR,10,'grid angle',IRET)
Call NCAPTC(CDFID,ANGID,'units',NCCHAR,7,'radians',IRET)
VRNG(1) = HMIN
VRNG(2) = HMAX
Write (*,*) 'hmin,hmax', HMIN, HMAX
Call NCAPT(CDFID,DEPTHID,'valid_range',NCFLOAT,2,VRNG,IRET)
*
Call NCAPTC(CDFID,TID,'long_name',NCCHAR,4,'Time',IRET)
Call NCAPTC(CDFID,TID,'units',NCCHAR,4,'days',IRET)
*
VRNG(1) = ELMIN
VRNG(2) = ELMAX
Write (*,*) 'elmin,elmax', ELMIN, ELMAX
Call NCAPT(CDFID,ELID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,ELID,'long_name',NCCHAR,9,'Elevation',IRET)
Call NCAPTC(CDFID,ELID,'units',NCCHAR,6,'meters',IRET)
FLOATVAL(1) = ELSC
Call NCAPT(CDFID,ELID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = ELOF
Call NCAPT(CDFID,ELID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
*
c If (PTHF.EQ.'y'.OR.PTHF.EQ.'Y') Then
* define heatflux
c DIMS(3) = TIMEDIM
c DIMS(2) = YPOSDIM
c DIMS(1) = XPOSDIM
c HFID = NCVDEF(CDFID,'heat_flux',NCSHORT,3,DIMS,IRET)
c VRNG(1) = HFMIN
c VRNG(2) = HFMAX
c Write (*,*) 'hfmin,hfmax', HFMIN, HFMAX
c Call NCAPT(CDFID,HFID,'valid_range',NCFLOAT,2,VRNG,IRET)
c Call NCAPTC(CDFID,HFID,'long_name',NCCHAR,9,'heat_flux',IRET)
c Call NCAPTC(CDFID,HFID,'units',NCCHAR,7,'watt/m2',IRET)
c FLOATVAL(1) = HFSC
c Call NCAPT(CDFID,HFID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
c FLOATVAL(1) = HFOF
c Call NCAPT(CDFID,HFID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
c End If
* define taux
DIMS(3) = TIMEDIM
DIMS(2) = YPOSDIM
DIMS(1) = XPOSDIM
TXID = NCVDEF(CDFID,'taux',NCSHORT,3,DIMS,IRET)
TYID = NCVDEF(CDFID,'tauy',NCSHORT,3,DIMS,IRET)
VRNG(1) = TAUMIN
VRNG(2) = TAUMAX
Write (*,*) 'TXmin,TXmax', TXMIN, TXMAX
Call NCAPTC(CDFID,TXID,'long_name',NCCHAR,16,'East Wind Stress',
* IRET)
Call NCAPTC(CDFID,TXID,'units',NCCHAR,10,'Newtons/m2',IRET)
Call NCAPTC(CDFID,TYID,'units',NCCHAR,10,'Newtons/m2',IRET)
Call NCAPTC(CDFID,TYID,'long_name',NCCHAR,17,
* 'North Wind Stress',IRET)
Call NCAPT(CDFID,TXID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPT(CDFID,TYID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPT(CDFID,TXID,'scale_factor',NCFLOAT,1,TAUSC,IRET)
Call NCAPT(CDFID,TYID,'scale_factor',NCFLOAT,1,TAUSC,IRET)
Call NCAPT(CDFID,TXID,'add_offset',NCFLOAT,1,TAUOF,IRET)
Call NCAPT(CDFID,TYID,'add_offset',NCFLOAT,1,TAUOF,IRET)
c PTCD='y'
c If (PTCD.EQ.'y'.OR.PTCD.EQ.'Y') Then
c define bottom drag
c DIMS(3) = TIMEDIM
c DIMS(2) = YPOSDIM
c DIMS(1) = XPOSDIM
c HFID = NCVDEF(CDFID,'cd',NCFLOAT,3,DIMS,IRET)
c VRNG(1) = 100.
c VRNG(2) = 0.
c Write (*,*) 'cdmin,cdmax', cdMIN, cdMAX
c Call NCAPT(CDFID,HFID,'valid_range',NCFLOAT,2,VRNG,IRET)
c Call NCAPTC(CDFID,HFID,'long_name',NCCHAR,24,
c & 'Bottom Drag Coefficient',IRET)
c Call NCAPTC(CDFID,HFID,'units',NCCHAR,4,'none',IRET)
c End If
* 4D Vars
DIMS(4) = TIMEDIM
DIMS(3) = ZPOSDIM
c
If (PTU.EQ.'Y'.OR.PTU.EQ.'y') Then
c u velocity
UID = NCVDEF(CDFID,'u',NCSHORT,4,DIMS,IRET)
VRNG(1) = VELMIN
VRNG(2) = VELMAX
Write (*,*) 'umin,umax', VELMIN, VELMAX, VRNG
Call NCAPT(CDFID,UID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,UID,'long_name',NCCHAR,11,'U1 Velocity',IRET
* )
Call NCAPTC(CDFID,UID,'units',NCCHAR,3,'m/s',IRET)
FLOATVAL(1) = VELSC
Call NCAPT(CDFID,UID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = VELOF
Call NCAPT(CDFID,UID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
End If
c
If (PTV.EQ.'Y'.OR.PTV.EQ.'y') Then
c v velocity
VID = NCVDEF(CDFID,'v',NCSHORT,4,DIMS,IRET)
VRNG(1) = VELMIN
VRNG(2) = VELMAX
Write (*,*) 'vmin,vmax', VELMIN, VELMAX, VRNG
Call NCAPT(CDFID,VID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,VID,'long_name',NCCHAR,11,'V1 Velocity',IRET
* )
Call NCAPTC(CDFID,VID,'units',NCCHAR,3,'m/s',IRET)
FLOATVAL(1) = VELSC
Call NCAPT(CDFID,VID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = VELOF
Call NCAPT(CDFID,VID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
End If
If (PTW.EQ.'Y'.OR.PTW.EQ.'y') Then
c w velocity
WID = NCVDEF(CDFID,'w',NCSHORT,4,DIMS,IRET)
VRNG(1) = WMIN
VRNG(2) = WMAX
Write (*,*) 'wmin,wmax', WMIN, WMAX, VRNG
Call NCAPT(CDFID,WID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,WID,'long_name',NCCHAR,10,'W velocity',IRET)
Call NCAPTC(CDFID,WID,'units',NCCHAR,3,'m/s',IRET)
FLOATVAL(1) = WSC
Call NCAPT(CDFID,WID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = WOF
Call NCAPT(CDFID,WID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
End If
If (PTC.EQ.'Y'.OR.PTC.EQ.'y') Then
c concentration
CCID = NCVDEF(CDFID,'conc',NCSHORT,4,DIMS,IRET)
VRNG(1) = CMIN
VRNG(2) = CMAX
Write (*,*) 'cmin,cmax', CMIN, CMAX, VRNG
Call NCAPT(CDFID,CCID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,CCID,'long_name',NCCHAR,13,'Concentration',
* IRET)
FLOATVAL(1) = CSC
Call NCAPT(CDFID,CCID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = COF
Call NCAPT(CDFID,CCID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
Call NCAPTC(CDFID,CCID,'units',NCCHAR,3,'ppt',IRET)
End If
If (PTS.EQ.'Y'.OR.PTS.EQ.'y') Then
c salinity
SALID = NCVDEF(CDFID,'salt',NCSHORT,4,DIMS,IRET)
VRNG(1) = SALMIN
VRNG(2) = SALMAX
Write (*,*) 'salmin,salmax', SALMIN, SALMAX, VRNG
Call NCAPT(CDFID,SALID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,SALID,'long_name',NCCHAR,8,'Salinity',IRET)
FLOATVAL(1) = SALSC
Call NCAPT(CDFID,SALID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = SALOF
Call NCAPT(CDFID,SALID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
Call NCAPTC(CDFID,SALID,'units',NCCHAR,3,'ppt',IRET)
End If
If (PTT.EQ.'Y'.OR.PTT.EQ.'y') Then
c temperature
TMPID = NCVDEF(CDFID,'temp',NCSHORT,4,DIMS,IRET)
VRNG(1) = TMPMIN
VRNG(2) = TMPMAX
Write (*,*) 'tmpmin,tmpmax', TMPMIN, TMPMAX, VRNG
Call NCAPT(CDFID,TMPID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,TMPID,'long_name',NCCHAR,11,'Temperature',
* IRET)
Call NCAPTC(CDFID,TMPID,'units',NCCHAR,7,'Celsius',IRET)
FLOATVAL(1) = TMPSC
Call NCAPT(CDFID,TMPID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = TMPOF
Call NCAPT(CDFID,TMPID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
End If
If (PTKM.EQ.'Y'.OR.PTKM.EQ.'y') Then
c vertical viscosity
KMID = NCVDEF(CDFID,'km',NCSHORT,4,DIMS,IRET)
VRNG(1) = KMMIN
VRNG(2) = KMMAX
Write (*,*) 'kmmin,kmmax', KMMIN, KMMAX, VRNG
Call NCAPT(CDFID,KMID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,KMID,'long_name',NCCHAR,18,
* 'vertical viscosity',IRET)
Call NCAPTC(CDFID,KMID,'units',NCCHAR,4,'m2/s',IRET)
FLOATVAL(1) = KMSC
Call NCAPT(CDFID,KMID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = KMOF
Call NCAPT(CDFID,KMID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
End If
If (PTAM.EQ.'Y'.OR.PTAM.EQ.'y') Then
c horizontal eddy viscosity
AMID = NCVDEF(CDFID,'am',NCSHORT,4,DIMS,IRET)
VRNG(1) = AMMIN
VRNG(2) = AMMAX
Write (*,*) 'ammin,ammax', AMMIN, AMMAX
Call NCAPT(CDFID,AMID,'valid_range',NCFLOAT,2,VRNG,IRET)
Call NCAPTC(CDFID,AMID,'long_name',NCCHAR,25,
* 'horizontal eddy viscosity',IRET)
Call NCAPTC(CDFID,AMID,'units',NCCHAR,4,'m2/s',IRET)
FLOATVAL(1) = AMSC
Call NCAPT(CDFID,AMID,'scale_factor',NCFLOAT,1,FLOATVAL,IRET)
FLOATVAL(1) = AMOF
Call NCAPT(CDFID,AMID,'add_offset',NCFLOAT,1,FLOATVAL,IRET)
End If
c
c Global attributes
Call NCAPTC(CDFID,NCGLOBAL,'experiment',NCCHAR,80,
* COM1,IRET)
c
c rps Hardwire start date to 1 Jan 1900 since
c POM doesn't know about Gregorian dates
IYR=1900
IMO=1
IDA=1
c rps
BASE_DATE(1) = IYR
BASE_DATE(2) = IMO
BASE_DATE(3) = IDA
Call NCAPT(CDFID,NCGLOBAL,'base_date',NCLONG,3,BASE_DATE,IRET)
c Call NCAPTC(CDFID,NCGLOBAL,'type_of_run',NCCHAR,10,TOR,IRET)
c Call NCAPT(CDFID,NCGLOBAL,'hormix_type',NCCHAR,10,HORZMIX,
c * IRET)
c Call NCAPT(CDFID,NCGLOBAL,'hormix_value',NCFLOAT,1,HORCON,
c * IRET)
c Call NCAPT(CDFID,NCGLOBAL,'vrtmix_type',NCCHAR,10,VERTMIX,
c * IRET)
c Call NCAPT(CDFID,NCGLOBAL,'vrtmix_value',NCFLOAT,1,UMOL,
c * IRET)
c Call NCAPT(CDFID,NCGLOBAL,'tlag',NCFLOAT,1,1./tlag,IRET)
c Call NCAPT(CDFID,NCGLOBAL,'tlax_w',NCFLOAT,1,TRLXW,IRET)
c Call NCAPT(CDFID,NCGLOBAL,'tlax_s',NCFLOAT,1,TRLXS,IRET)
c Call NCAPT(CDFID,NCGLOBAL,'tlax_e',NCFLOAT,1,TRLXE,IRET)
c Call NCAPT(CDFID,NCGLOBAL,'tlax_n',NCFLOAT,1,TRLXN,IRET)
Call NCENDF(CDFID,IRET)
c=========================================================================
c END OF DEFINE MODE
c=========================================================================
c store easting
Write (*,*) 'storing easting'
CORNER(1) = 1
CORNER(2) = 1
EDGES(1) = IM
EDGES(2) = JM
Call NCVPT(CDFID,XID,CORNER,EDGES,XGRID,IRET)
c store northing
Write (*,*) 'storing northing'
Call NCVPT(CDFID,YID,CORNER,EDGES,YGRID,IRET)
c store h1
Call NCVPT(CDFID,H1ID,CORNER,EDGES,DX,IRET)
Write (*,*) 'storing dx'
c store h2
Call NCVPT(CDFID,H2ID,CORNER,EDGES,DY,IRET)
Write (*,*) 'storing dy'
c store depth
do j=1,jm
do i=1,im
if(fsm(i,j).eq.0.0)then
c ECOMVIZ determines land cells by looking for depth=-99999.
htmp(i,j)=-99999.
else
htmp(i,j)=h(i,j)
endif
enddo
enddo
Write (*,*) 'storing depth'
Call NCVPT(CDFID,DEPTHID,CORNER,EDGES,HTMP,IRET)
c store grid angles
Write (*,*) 'storing grid angles'
Call NCVPT(CDFID,ANGID,CORNER,EDGES,ANG,IRET)
c store stretched vertical coordinate levels (sigma-levels)
Write (*,*) 'storing sigma levels'
CORNER(1) = 1
EDGES(1) = KB
Call NCVPT(CDFID,SIGMAID,CORNER,EDGES,Z,IRET)
Else
CDFID = NCOPN('pom97.cdf',NCWRITE,IRET)
IFILL = NCSFIL(CDFID,NCNOFILL,IRET)
End If
ICDF = ICDF + 1
c store time
CORNER(1) = ICDF
Write (*,*) 'storing time'
TID = NCVID(CDFID,'time',IRET)
Call NCVPT1(CDFID,TID,CORNER,TMIDDLE,IRET)
c store elevation
CORNER(1) = 1
CORNER(2) = 1
CORNER(3) = ICDF
EDGES(1) = IM
EDGES(2) = JM
EDGES(3) = 1
Do J = 1, JM
Do I = 1, IM
WS0(I,J) = NINT((ELB(I,J)-ELOF)/ELSC)
End Do
End Do
Write (*,*) 'storing elevations'
ELID = NCVID(CDFID,'elev',IRET)
Call NCVPT(CDFID,ELID,CORNER,EDGES,WS0,IRET)
c store wind stress components
Write (*,*) 'storing wind stress'
TXID = NCVID(CDFID,'taux',IRET)
Do J = 1, JM
Do I = 1, IM
WS0(I,J) = NINT((-WUSURF(I,J)*1000.-TAUOF)/TAUSC)
End Do
End Do
TYID = NCVID(CDFID,'tauy',IRET)
Call NCVPT(CDFID,TXID,CORNER,EDGES,WS0,IRET)
Do J = 1, JM
Do I = 1, IM
WS0(I,J) = NINT((-WVSURF(I,J)*1000.-TAUOF)/TAUSC)
End Do
End Do
Call NCVPT(CDFID,TYID,CORNER,EDGES,WS0,IRET)
c
c If (PTCD.EQ.'y'.OR.PTCD.EQ.'Y') Then
c CDID = NCVID(CDFID,'cd',IRET)
c Call NCVPT(CDFID,cdID,CORNER,EDGES,ARCCD,IRET)
c endif
CORNER(3) = 1
CORNER(4) = ICDF
EDGES(3) = KB
EDGES(4) = 1
c store concentration
c If (PTC.EQ.'Y'.OR.PTC.EQ.'y') Then
c Do I = 1, IM
c Do J = 1, JM
c Do K = 1, KB
c WS1(I,J,K) = NINT((MAX(ARCC(I,J,K),0.)-COF)/CSC)
c End Do
c End Do
c End Do
c Write (*,*) 'storing concentration'
c CCID = NCVID(CDFID,'conc',IRET)
c Call NCVPT(CDFID,CCID,CORNER,EDGES,WS1,IRET)
c End If
c store heatflux
c If (PTHF.EQ.'y'.OR.PTHF.EQ.'Y') Then
c CORNER(1) = 1
c CORNER(2) = 1
c CORNER(3) = ICDF
c EDGES(1) = IM
c EDGES(2) = JM
c EDGES(3) = 1
c Do I = 1, IM
c Do J = 1, JM
c WS0(I,J) = NINT((ARCHF(I,J)*(-4.186E6)-HFOF)/HFSC)
c End Do
c End Do
c Write (*,*) 'storing heat flux'
c HFID = NCVID(CDFID,'heat_flux',IRET)
c Call NCVPT(CDFID,HFID,CORNER,EDGES,WS0,IRET)
c End If
CORNER(3) = 1
CORNER(4) = ICDF
EDGES(3) = KB
EDGES(4) = 1
c store salinity
If (PTS.EQ.'Y'.OR.PTS.EQ.'y') Then
Do I = 1, IM
Do J = 1, JM
Do K = 1, KB
WS1(I,J,K) = NINT((SB(I,J,K)+SBIAS-SALOF)/SALSC)
End Do
End Do
End Do
Write (*,*) 'storing salinity'
SALID = NCVID(CDFID,'salt',IRET)
Call NCVPT(CDFID,SALID,CORNER,EDGES,WS1,IRET)
End If
c store temperature
If (PTT.EQ.'Y'.OR.PTT.EQ.'y') Then
Do I = 1, IM
Do J = 1, JM
Do K = 1, KB
WS1(I,J,K) = NINT((TB(I,J,K)+TBIAS-TMPOF)/TMPSC)
End Do
End Do
End Do
Write (*,*) 'storing temperature'
TMPID = NCVID(CDFID,'temp',IRET)
Call NCVPT(CDFID,TMPID,CORNER,EDGES,WS1,IRET)
End If
c store vertical viscosity
If (PTKM.EQ.'Y'.OR.PTKM.EQ.'y') Then
Do K = 1, KB
Do J = 1, JM
Do I = 1, IM
WS1(I,J,K) = NINT((KM(I,J,K)-KMOF)/KMSC)
End Do
End Do
End Do
Write (*,*) 'storing vertical viscosity'
KMID = NCVID(CDFID,'km',IRET)
Call NCVPT(CDFID,KMID,CORNER,EDGES,WS1,IRET)
End If
c store horizontal viscosity
If (PTAM.EQ.'Y'.OR.PTAM.EQ.'y') Then
Do I = 1, IM
Do J = 1, JM
Do K = 1, KB
WS1(I,J,K) = NINT((AAM(I,J,K)-AMOF)/AMSC)
End Do
End Do
End Do
Write (*,*) 'storing horizontal viscosity'
AMID = NCVID(CDFID,'am',IRET)
Call NCVPT(CDFID,AMID,CORNER,EDGES,WS1,IRET)
End If
c store velocity
If (PTU.EQ.'Y'.OR.PTU.EQ.'y') Then
Do I = 1, IM
Do J = 1, JM
Do K = 1, KB
WS1(I,J,K) = NINT((U(I,J,K)-VELOF)/VELSC)
End Do
End Do
End Do
Write (*,*) 'storing U component of velocity'
UID = NCVID(CDFID,'u',IRET)
Call NCVPT(CDFID,UID,CORNER,EDGES,WS1,IRET)
End If
If (PTV.EQ.'Y'.OR.PTV.EQ.'y') Then
Do I = 1, IM
Do J = 1, JM
Do K = 1, KB
WS1(I,J,K) = NINT((V(I,J,K)-VELOF)/VELSC)
End Do
End Do
End Do
Write (*,*) 'storing V component of velocity'
VID = NCVID(CDFID,'v',IRET)
Call NCVPT(CDFID,VID,CORNER,EDGES,WS1,IRET)
End If
If (PTW.EQ.'Y'.OR.PTW.EQ.'y') Then
Do I = 1, IM
Do J = 1, JM
Do K = 1, KB
WS1(I,J,K) = NINT((W(I,J,K)-WOF)/WSC)
End Do
End Do
End Do
Write (*,*) 'storing W component of velocity'
WID = NCVID(CDFID,'w',IRET)
Call NCVPT(CDFID,WID,CORNER,EDGES,WS1,IRET)
End If
Call NCCLOS(CDFID,IRET)
End
]]>