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