c c --- profout first (before PROGRAM) to fix an AIX xlf95 bug. c subroutine profout(array,zz, plonj,platj,distj,jn, & artype,yrflag,time3,nstep,iexpt, & name,namec,names,units, kz, ncfile) use netcdf ! NetCDF Interface implicit none c character*(*) name,namec,names,units,ncfile integer jn,artype,yrflag,nstep,iexpt,kz double precision time3(3) real array(kz,jn),zz(kz), & plonj(jn),platj(jn),distj(jn) c c write out a 2-d z-level array to a netCDF file. c c 2-d array size must be identical in all calls. c artype,yrflag,time3 must be identical in all calls. c c the NetCDF filename is taken from ncfile. c the NetCDF title, institution and section_name (if any) are taken c from environment variables CDF_TITLE, CDF_INST and CDF_SECTION. c If environment variables MERSEA_B_DATE and MERSEA_B_TYPE exist, c then the output format is in MERSEA format (i.e. no time axis). c c This version needs version 3.5 of the NetCDF library, from: c http://www.unidata.ucar.edu/packages/netcdf/ c real :: array_t(jn,kz) !transposed array real :: array_0(jn) !surface array c integer :: ncfileID, status, varID, latVarID,lonVarID integer :: StaDimID,StaVarID, & DepDimID,DepVarID integer :: MTDimID,MTVarID,datVarID character :: ncenv*240 character :: Ename*6 c logical :: lopen,lexist,lmersea,ltrack,lstation,llon,llat integer :: i,j,k,l,iyear,month,iday,ihour, & iyrms,monms,idms,ihrms real :: hmin(999),hmax(999) double precision :: dt,dt0,year c character*81, save :: label = ' ' integer, save :: ioinit = -1 double precision, save :: date = 0.d0 double precision, save :: cell = 0.d0 real, parameter :: fill_value = 2.0**100 c character cmonth(12)*3 data cmonth/'jan','feb','mar','apr','may','jun', & 'jul','aug','sep','oct','nov','dec'/ c ncenv = ' ' call getenv('MERSEA_B_TYPE',ncenv) lmersea = ncenv.ne.' ' lstation = distj(jn).eq.0.0 ltrack = distj(jn).gt.0.0 llon = distj(jn).lt.0.0 .and. distj(jn).gt.-1.5 llat = distj(jn).lt.0.0 .and. distj(jn).le.-1.5 * write(6,*) * write(6,*) 'MERSEA_B_TYPE: ',trim(ncenv) * write(6,*) 'lmersea = ',lmersea * write(6,*) 'lstation = ',lstation * write(6,*) 'ltrack = ',ltrack * write(6,*) 'llon = ',llon * write(6,*) 'llat = ',llat * write(6,*) 'jn,kz = ',jn,kz * write(6,*) c if (ioinit.eq.-1) then ioinit = 0 c c initialize label. c if (yrflag.eq.0) then year = 360.0d0 elseif (yrflag.lt.3) then year = 366.0d0 else year = 365.25d0 endif call fordate(time3(3),yrflag, iyear,month,iday,ihour) date = (iday + 100 * month + 10000 * iyear) + & (time3(3)-int(time3(3))) if (artype.eq.1) then write (label(51:72),123) cmonth(month),iday,iyear,ihour elseif (artype.eq.2 .and. time3(2)-time3(1).lt.1.1) then !daily mean write (label(51:72),223) cmonth(month),iday,iyear else ! mean or sdev archive write(6,*) 'time3 = ',time3 dt = 0.5*(time3(2)-time3(1))/(nstep-1) if (yrflag.eq.0) then dt0 = 15.0 elseif (yrflag.eq.1) then dt0 = 15.25 elseif (yrflag.eq.2) then dt0 = 0.0 else dt0 = 0.0 endif cell = (time3(2)+dt+dt0) - (time3(1)-dt+dt0) if (artype.eq.2) then write(label(51:72),114) ' mean: ',(time3(1)-dt+dt0)/year, & (time3(2)+dt+dt0)/year else write(label(51:72),114) ' sdev: ',(time3(1)-dt+dt0)/year, & (time3(2)+dt+dt0)/year endif endif write (label(73:81),115) iexpt/10,mod(iexpt,10),'H' 123 format (' ',a3,i3.2,',',i5.4,i3.2,'Z ') 223 format (' ',a3,i3.2,',',i5.4,' MEAN ') 114 format (a7,f7.2,'-',f7.2) 115 format (' [',i2.2,'.',i1.1,a1,']') endif !initialization c c NetCDF I/O c call ncrange(array, 1,jn,kz, fill_value, hmin(1),hmax(1)) c inquire(file= ncfile, exist=lexist) if (.not.lexist) then ! create a new NetCDF and write data to it ! netcdf-4 classic model, netcdf version 4.3 and later call nchek('nf90_create', & nf90_create(trim(ncfile), & or(nf90_clobber, & or(nf90_hdf5, & nf90_classic_model)), & ncfileID)) ! define the dimensions if (.not.lmersea) then call nchek("nf90_def_dim", & nf90_def_dim(ncfileID, & "MT", nf90_unlimited,MTDimID)) endif if (llon) then call nchek("nf90_def_dim", & nf90_def_dim(ncfileID,"longitude",jn,StaDimID)) elseif (llat) then call nchek("nf90_def_dim", & nf90_def_dim(ncfileID,"latitude", jn,StaDimID)) elseif (ltrack) then call nchek("nf90_def_dim", & nf90_def_dim(ncfileID,"track", jn,StaDimID)) else !station call nchek("nf90_def_dim", & nf90_def_dim(ncfileID,"station", jn,StaDimID)) endif call nchek("nf90_def_dim", & nf90_def_dim(ncfileID,"depth", kz,DepDimID)) ! create the global attributes call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "Conventions", & "CF-1.0")) ncenv = ' ' call getenv('CDF_TITLE',ncenv) if (ncenv.eq.' ') then ncenv = "HYCOM" endif call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "title", & trim(ncenv))) ncenv = ' ' call getenv('CDF_INST',ncenv) if (ncenv.ne.' ') then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "institution", & trim(ncenv))) endif ncenv = ' ' call getenv('CDF_SECTION',ncenv) if (ncenv.ne.' ') then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "section_name", & trim(ncenv))) endif if (artype.eq.1) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "source", & "HYCOM archive file")) elseif (artype.eq.2) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "source", & "HYCOM mean archive file")) else call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "source", & "HYCOM std. archive file")) endif call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "experiment", & label(75:78))) call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "history", & "hycom_profile2z_nc")) if (lmersea) then if (artype.eq.2) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "field_type", & "daily average")) else call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "field_type", & "instantaneous")) endif write(ncenv, & '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":00:00")') & iyear,month,iday,ihour call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "field_date", & trim(ncenv))) ncenv = ' ' call getenv('MERSEA_B_DATE',ncenv) if (ncenv.eq.'TODAY') then write(ncenv, & '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":00:00")') & iyear,month,iday,ihour endif if (ncenv.ne.' ') then read(ncenv,'(i4,1x,i2,1x,i2,1x,i2)') & iyrms,monms,idms,ihrms if (iyrms.lt.iyear) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "forecast")) elseif (iyrms.gt.iyear) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "hindcast")) else !iyrms.eq.iyear if (monms.lt.month) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "forecast")) elseif (monms.gt.month) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "hindcast")) else !monms.eq.month if (idms.lt.iday) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "forecast")) elseif (idms.gt.iday) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "hindcast")) else !idms.eq.iday if (ihrms.lt.ihour) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "forecast")) elseif (ihrms.gt.ihour) then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "hindcast")) else !ihrms.eq.ihour call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "forecast_type", & "nowcast")) endif !ihrms endif !idms endif !monms endif !iyrms call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "bulletin_date", & trim(ncenv))) endif ncenv = ' ' call getenv('MERSEA_B_TYPE',ncenv) if (ncenv.ne.' ') then call nchek("nf90_put_att (global)", & nf90_put_att(ncfileID,nf90_global, & "bulletin_type", & trim(ncenv))) endif endif !lmersea ! create the variables and attributes if (.not.lmersea) then call nchek("nf90_def_var (MT)", & nf90_def_var(ncfileID,"MT", nf90_double, & MTDimID,MTVarID)) if (yrflag.eq.0) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "long_name", & "model time")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "units", & "days since 0001-01-01 00:00:00")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "calendar", & "360_day")) elseif (yrflag.eq.1) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "long_name", & "model time")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "units", & "days since 0001-01-16 00:00:00")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "calendar", & "366_day")) elseif (yrflag.eq.2) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "long_name", & "model time")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "units", & "days since 0001-01-01 00:00:00")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "calendar", & "366_day")) else call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "long_name", & "time")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "units", & "days since 1900-12-31 00:00:00")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "calendar", & "standard")) endif call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "axis","T")) call nchek("nf90_def_var (Date)", & nf90_def_var(ncfileID,"Date", nf90_double, & MTDimID,datVarID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "long_name", & "date")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "units", & "day as %Y%m%d.%f")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "C_format", & "%13.4f")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "FORTRAN_format", & "(f13.4)")) if (artype.eq.2) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "cell_methods", & "mean")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "cell_extent", & cell)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "cell_methods", & "mean")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "cell_extent", & cell)) elseif (artype.eq.3) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "cell_methods", & "standard_deviation")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,MTVarID, & "cell_extent", & cell)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "cell_methods", & "standard_deviation")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,datVarID, & "cell_extent", & cell)) endif endif !.not.lmersea call nchek("nf90_def_var (depth)", & nf90_def_var(ncfileID,"depth", nf90_float, & DepDimID,DepVarID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,DepVarID, & "units","m")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,DepVarID, & "positive","down")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,DepVarID, & "axis","Z")) if (ltrack) then call nchek("nf90_def_var (track)", & nf90_def_var(ncfileID,"track", nf90_float, & StaDimID,StaVarID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,StaVarID, & "units","km")) if (jn.ne.1 .and. & abs( (distj(jn)-distj(1))- & (distj( 2)-distj(1))*(jn-1) ).lt.0.01) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,StaVarID, & "point_spacing","even")) !ferret endif elseif (lstation) then call nchek("nf90_def_var (station)", & nf90_def_var(ncfileID,"station", nf90_int, & StaDimID,StaVarID)) endif ! longitude and latitude call nchek("nf90_def_var (latitude)", & nf90_def_var(ncfileID,"latitude", nf90_float, & StaDimID,latVarID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,latVarID, & "standard_name","latitude")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,latVarID, & "units","degrees_north")) if (llat) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,latVarID, & "axis","Y")) endif call nchek("nf90_put_att", & nf90_put_att(ncfileID,latVarID, & "C_format","%8.3f")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,latVarID, & "FORTRAN_format","(f8.3)")) call nchek("nf90_def_var (longitude)", & nf90_def_var(ncfileID,"longitude", nf90_float, & StaDimID,lonVarID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,lonVarID, & "standard_name","longitude")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,lonVarID, & "units","degrees_east")) if (llon) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,lonVarID, & "axis","X")) endif call nchek("nf90_put_att", & nf90_put_att(ncfileID,lonVarID, & "C_format","%8.3f")) call nchek("nf90_put_att", & nf90_put_att(ncfileID,lonVarID, & "FORTRAN_format","(f8.3)")) ! model 2Z variable if (kz.gt.1) then if (lmersea) then * write(6,*) 'StaDimID = ',StaDimID * write(6,*) 'DepDimID = ',DepDimID call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & (/StaDimID, DepDimID/), & varID)) if (llon) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "latitude")) elseif (llat) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude")) else !track/station call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude")) endif !lon:lat:else else !.not.lmersea call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & (/StaDimID, & DepDimID, MTDimID/), & varID)) if (llon) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "latitude Date")) elseif (llat) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude Date")) else !track/station call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude Date")) endif !lon:lat:else endif !mersea:else else !kz.eq.1 if (lmersea) then * write(6,*) 'StaDimID = ',StaDimID call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & StaDimID, varID)) if (llon) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "latitude")) elseif (llat) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude")) else !track/station call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude")) endif !lon:lat:else else !.not.lmersea call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & (/StaDimID, MTDimID/), & varID)) if (llon) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "latitude Date")) elseif (llat) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude Date")) else !track/station call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude Date")) endif !lon:lat:else endif !mersea:else endif !kz.gt.1:else if (names.ne." ") then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "standard_name",trim(names))) endif call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID,"units",trim(units))) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "_FillValue",fill_value)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "valid_range", & (/hmin(1), hmax(1)/))) if (artype.eq.1) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "long_name", & trim(name)//label(73:81))) else call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "long_name", & trim(name)//label(51:55)//label(73:81))) endif ! leave def mode call nchek("nf90_enddef", & nf90_enddef(ncfileID)) ! write data into coordinate variables if (.not.lmersea) then call nchek("nf90_put_var", & nf90_put_var(ncfileID,MTVarID, time3(3))) call nchek("nf90_put_var", & nf90_put_var(ncfileID,datVarID,date )) endif if (ltrack) then call nchek("nf90_put_var", & nf90_put_var(ncfileID,StaVarID,distj(:))) elseif (lstation) then call nchek("nf90_put_var", & nf90_put_var(ncfileID,StaVarID, & (/ (j,j=1,jn) /))) endif call nchek("nf90_put_var", & nf90_put_var(ncfileID,DepVarID,zz(:) )) call nchek("nf90_put_var", & nf90_put_var(ncfileID,latVarID,platj(:))) call nchek("nf90_put_var", & nf90_put_var(ncfileID,lonVarID,plonj(:))) ! write to model variable if (kz.gt.1) then do j= 1,jn do k= 1,kz array_t(j,k) = array(k,j) enddo enddo call nchek("nf90_put_var", & nf90_put_var(ncfileID,varID,array_t(:,:))) else do j= 1,jn array_0(j) = array(1,j) enddo call nchek("nf90_put_var", & nf90_put_var(ncfileID,varID,array_0(:))) endif ! close NetCDF file call nchek("nf90_close", & nf90_close(ncfileID)) else c c Write data to existing NetCDF file c ! open NetCDF file call nchek("nf90_open", & nf90_open(trim(ncfile), nf90_write, ncfileID)) ! get dimension ID's if (.not.lmersea) then call nchek("nf90_inq_dimid", & nf90_inq_dimid(ncfileID,"MT",MTDimID)) endif call nchek("nf90_inq_dimid", & nf90_inq_dimid(ncfileID,"depth", DepDimID)) if (llon) then call nchek("nf90_inq_dimid", & nf90_inq_dimid(ncfileID,"longitude",StaDimID)) elseif (llat) then call nchek("nf90_inq_dimid", & nf90_inq_dimid(ncfileID,"latitude", StaDimID)) elseif (ltrack) then call nchek("nf90_inq_dimid", & nf90_inq_dimid(ncfileID,"track", StaDimID)) elseif (lstation) then call nchek("nf90_inq_dimid", & nf90_inq_dimid(ncfileID,"station", StaDimID)) endif ! switch to define mode call nchek("nf90_redef", & nf90_redef(ncfileID)) ! define new variable if (kz.gt.1) then if (lmersea) then call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & (/StaDimID, DepDimID/), & varID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude")) else !.not.lmersea call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & (/StaDimID, DepDimID, MTDimID/), & varID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude Date")) endif !lmersea:else else !kz.eq.1 if (lmersea) then call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & StaDimID, & varID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude")) else !.not.lmersea call nchek("nf90_def_var (namec)", & nf90_def_var(ncfileID,trim(namec),nf90_float, & (/StaDimID, MTDimID/), & varID)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "coordinates", & "longitude latitude Date")) endif !lmersea:else endif !kz if (names.ne." ") then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "standard_name",trim(names))) endif call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID,"units",trim(units))) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "_FillValue",fill_value)) call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "valid_range", & (/hmin(1), hmax(1)/))) if (artype.eq.1) then call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "long_name", & trim(name)//label(73:81))) else call nchek("nf90_put_att", & nf90_put_att(ncfileID,varID, & "long_name", & trim(name)//label(51:55)//label(73:81))) endif ! leave define mode call nchek("nf90_enddef", & nf90_enddef(ncfileID)) ! inquire variable ID call nchek("nf90_inq_varid", & nf90_inq_varid(ncfileID,trim(namec),varID)) !write values into array if (kz.gt.1) then do j= 1,jn do k= 1,kz array_t(j,k) = array(k,j) enddo enddo call nchek("nf90_put_var", & nf90_put_var(ncfileID,varID,array_t(:,:))) else do j= 1,jn array_0(j) = array(1,j) enddo call nchek("nf90_put_var", & nf90_put_var(ncfileID,varID,array_0(:))) endif ! close file call nchek("nf90_close", & nf90_close(ncfileID)) endif write(6,'(a49,a,2g15.6)') & trim(name)//label(51:81),':',hmin(1),hmax(1) * call flush(6) return contains subroutine nchek(cnf90,status) implicit none c character*(*), intent(in) :: cnf90 integer, intent(in) :: status c c subroutine to handle NetCDF errors c * write(6,'(a)') trim(cnf90) * call flush(6) * if (status /= nf90_noerr) then write(6,'(/a)') 'error in profout - from NetCDF library' write(6,'(a/)') trim(cnf90) write(6,'(a/)') trim(nf90_strerror(status)) * call flush(6) stop end if end subroutine nchek end subroutine profout subroutine ncrange(h,ii,jj,kk, fill_value, hmin,hmax) implicit none c integer, intent(in ) :: ii,jj,kk real, intent(in ) :: h(ii,jj,kk),fill_value real, intent(out) :: hmin,hmax c c return range of array, ignoring fill_value c integer i,j,k real hhmin,hhmax c hhmin = abs(fill_value) hhmax = -abs(fill_value) do k= 1,kk do j= 1,jj do i= 1,ii if (h(i,j,k).ne.fill_value) then hhmin = min(hhmin,h(i,j,k)) hhmax = max(hhmax,h(i,j,k)) endif enddo enddo enddo hmin = hhmin hmax = hhmax end subroutine ncrange PROGRAM PROFILE IMPLICIT NONE C C hycom_profile2z_nc - Usage: hycom_profile2z_nc archv.a list.txt z.txt archz.nc [itype [bot]] C C produce a netCDF z-level profile file at a list of points C C archv.a is assumed to be an HYCOM archive data file, with companion C header file archv.b. C list.txt contains a list of points, one per line C z.txt contains a list of z-depths, one per line C archz.nc will be the output z-level NetCDF profile file C itype is the interpolation type (default 1) C =-2; piecewise quadratic cross cell C =-1; piecewise linear across cell C = 0; piecewise constant across cell C = 1; linear interpolation between cell centers C bot ignore layers within bot of the bottom (default 0.0) C C Each line of list.txt should contain three values: x,y,dist. C The x,y values are the p-grid location location of the point, C note that hycom/bin/hycom_lonlat2xy will convert lon,lat to x,y. C The dist value is the distance from the first point in km, C or 0.0 if the points do not represent a track or transect. C To use lon (lat) as the standard axis set dist to -1.0 (-2.0). C C The files regional.grid.[ab] for this domain must be in the current C directory. They are used to calculate lon,lat values for each point. C C The following environment variables control the netCDF output: C C CDF_TITLE - NetCDF title C CDF_INST - NetCDF institution C CDF_SECTION - NetCDF section_name C MERSEA_B_DATE - MERSEA Bulletin date C MERSEA_B_TYPE - MERSEA Bulletin type C C If environment variables MERSEA_B_DATE and MERSEA_B_TYPE exist, C then the output format is in MERSEA format (i.e. no time axis). C C Primarily for MERSEA, so only T,S,u,v (and fsd,sst,sss) are output. C C this version for "serial" Unix systems. C C Alan J. Wallcraft, Naval Research Laboratory, July 2003. C REAL*4 ONEM,SPVAL PARAMETER (ONEM=9806.0, SPVAL=2.0**100) C REAL*4, ALLOCATABLE :: A(:,:),AA(:,:) REAL*4, ALLOCATABLE :: XP(:),YP(:),DP(:) REAL*4, ALLOCATABLE :: PLON(:),PLAT(:),PANG(:) REAL*4, ALLOCATABLE :: P(:,:),TK(:,:),SK(:,:), + UK(:,:),VK(:,:), + UB(:),VB(:),FSD(:),SST(:),SSS(:) REAL*4, ALLOCATABLE :: ZZ(:),TZ(:,:),SZ(:,:), + UZ(:,:),VZ(:,:) REAL*4 :: PAD(4096) C INTEGER IARGC INTEGER NARG CHARACTER*240 CARG C LOGICAL LSTERIC,LSEAICE,LFATAL INTEGER IDM,JDM,KDM,KTR,NSURF,NLAY,NDIF,NTRC, + ARTYPE,IEXPT,YRFLAG,SIGVER,I_ARCH(17) INTEGER ITYPE,IDEBUG,KZ,NPAD,N,NX REAL*4 DUMZ,FLAG,THBASE,BOT INTEGER NSTEP REAL*8 TIME3(3) CHARACTER*240 CFILEA,CFILEB,CFILEL,CFILEP,CFILEZ,CFILEVS,CFORMAT C CHARACTER*18 CASN LOGICAL LDONE INTEGER I,II,IP,IPP1,IPP2,J,JJ,JP,JPP1,JPP2,K,LANDF,LANDL, + KREC,KREC0,IOS,NRECL REAL*4 A00,A01,A10,A11,AXY,DX,DY #ifdef CRAY INTEGER*8 IU8,IOS8 #endif C C READ ARGUMENTS. C NARG = IARGC() C IF (NARG.EQ.4) THEN CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEL) CALL GETARG(3,CFILEZ) CALL GETARG(4,CFILEP) ITYPE = 1 IDEBUG = 0 BOT = 0.0 ELSEIF (NARG.EQ.5) THEN CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEL) CALL GETARG(3,CFILEZ) CALL GETARG(4,CFILEP) CALL GETARG(5,CARG) READ(CARG,*) ITYPE IDEBUG = 0 BOT = 0.0 ELSEIF (NARG.EQ.6) THEN CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEL) CALL GETARG(3,CFILEZ) CALL GETARG(4,CFILEP) CALL GETARG(5,CARG) READ(CARG,*) ITYPE CALL GETARG(6,CARG) READ(CARG,*) BOT IDEBUG = 0 ELSEIF (NARG.EQ.7) THEN !undocumented debug option CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEL) CALL GETARG(3,CFILEZ) CALL GETARG(4,CFILEP) CALL GETARG(5,CARG) READ(CARG,*) ITYPE CALL GETARG(6,CARG) READ(CARG,*) BOT CALL GETARG(7,CARG) READ(CARG,*) IDEBUG !>0 print this station ELSE WRITE(6,"(2a)") + 'Usage: hycom_profile2z_nc', + ' archv.a list.txt z.txt archz.nc [itype [bot]]' CALL EXIT(1) STOP ENDIF C C EXTRACT MODEL PARAMETERS FROM ".b" FILE. C cfilevs = ' ' call getenv('ARCHVS',cfilevs) if (cfilevs.eq.' ') then CFILEB = CFILEA(1:LEN_TRIM(CFILEA)-1) // 'b' CALL READ_BS(CFILEB, + IEXPT,YRFLAG,IDM,JDM,KDM,NSURF,NLAY,NDIF,NTRC, + LSTERIC,ARTYPE,SIGVER,THBASE,TIME3,NSTEP) I_ARCH( :) = -1 !turn off i_arch else CFILEB = CFILEA(1:LEN_TRIM(CFILEA)-1) // 'b' CALL READ_BSCS(CFILEB,CFILEVS,I_ARCH, + IEXPT,YRFLAG,IDM,JDM,KDM,NSURF,NLAY,NDIF,NTRC, + LSEAICE,LSTERIC,ARTYPE,SIGVER,THBASE,TIME3,NSTEP) endif * write(6,*) 'IDM = ',idm * write(6,*) 'JDM = ',jdm * write(6,*) 'KDM = ',kdm C CALL SIG_I(SIGVER) C C FIND THE NUMBER OF PROFILES AND THE NUMBER OF Z-LEVELS. C OPEN(UNIT=31, FILE=CFILEL, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(CFILEL) WRITE(6,*) 'ios = ',ios CALL EXIT(3) STOP ENDIF C DO N= 1,99999 READ(31,*,IOSTAT=IOS) DX,DY,DUMZ IF (IOS.NE.0) THEN EXIT ENDIF ENDDO !n N = N-1 REWIND(31) * write(6,*) 'N = ',n C OPEN(UNIT=32, FILE=CFILEZ, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(CFILEZ) WRITE(6,*) 'ios = ',ios CALL EXIT(3) STOP ENDIF C DO KZ= 1,99999 READ(32,*,IOSTAT=IOS) DUMZ IF (IOS.NE.0) THEN EXIT ENDIF ENDDO !kz KZ = KZ-1 REWIND(32) * write(6,*) 'KZ = ',kz C C ALLOCATE ARRAYS. C NPAD = 4096 - MOD(IDM*JDM,4096) IF (NPAD.EQ.4096) THEN NPAD = 0 ENDIF C ALLOCATE( A(IDM,JDM), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ', + IDM*JDM,' words' CALL EXIT(2) STOP ENDIF ALLOCATE( AA(IDM+1,JDM), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ', + (IDM+1)*JDM,' words' CALL EXIT(2) STOP ENDIF ALLOCATE( P(0:KDM,N), + TK(1:KDM,N), + SK(1:KDM,N), + UK(1:KDM,N), + VK(1:KDM,N), + ZZ(1:KZ ), + TZ(1:KZ, N), + SZ(1:KZ, N), + UZ(1:KZ, N), + VZ(1:KZ, N), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ', + '(1:kz,1:n) arrays' CALL EXIT(2) STOP ENDIF ALLOCATE( UB( N), + VB( N), + FSD( N), + SST( N), + SSS( N), + PLON( N), + PLAT( N), + PANG( N), + XP( N), + YP( N), + DP( N), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ', + '(1:n) arrays' CALL EXIT(2) STOP ENDIF C C READ THE Z FILE. C DO K= 1,KZ READ(32,*) ZZ(K) IF (IOS.NE.0) THEN ENDIF ENDDO !k CLOSE(32) C IF (KZ.EQ.1 .AND. ZZ(1).EQ.0.0) THEN C C SURFACE ONLY. C KDM = 1 ENDIF C C READ PROFILE LOCATIONS C DO I= 1,N READ(31,*,IOSTAT=IOS) XP(I),YP(I),DP(I) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t read ',TRIM(CFILEL) WRITE(6,*) 'ios = ',ios WRITE(6,*) 'i,n = ',i,n CALL EXIT(3) STOP ENDIF ENDDO !i CLOSE(31) C C OPEN "regional.grid.a" FILE. C IF (NPAD.EQ.0) THEN INQUIRE( IOLENGTH=NRECL) A ELSE INQUIRE( IOLENGTH=NRECL) A,PAD(1:NPAD) ENDIF * write(6,*) 'nrecl = ',nrecl #ifdef CRAY #ifdef t3e IF (MOD(NRECL,4096).EQ.0) THEN WRITE(CASN,8000) NRECL/4096 8000 FORMAT('-F cachea:',I4.4,':1:0') IU8 = 11 CALL ASNUNIT(IU8,CASN,IOS8) IF (IOS8.NE.0) THEN WRITE(6,*) 'Error: can''t asnunit 11' WRITE(6,*) 'ios = ',ios8 WRITE(6,*) 'casn = ',casn CALL EXIT(5) STOP ENDIF ENDIF #else CALL ASNUNIT(11,'-F syscall -N ieee',IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t asnunit 11' WRITE(6,*) 'ios = ',ios CALL EXIT(5) STOP ENDIF #endif #endif OPEN(UNIT=11, FILE="regional.grid.a", + FORM='UNFORMATTED', STATUS='OLD', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open regional.grid.a' WRITE(6,*) 'ios = ',ios WRITE(6,*) 'nrecl = ',nrecl CALL EXIT(3) STOP ENDIF C READ(11,REC=1,IOSTAT=IOS) A !plon #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read regional.grid.a' WRITE(6,*) 'ios = ',ios CALL EXIT(5) STOP ENDIF AA(1:IDM,:) = A(:,:) DO J= 1,JDM A00 = MOD(AA(1,J) - AA(IDM,J),360.0) IF (A00.LT.-180.0) THEN A00 = A00 + 360.0 ELSEIF (A00.GT. 180.0) THEN A00 = A00 - 360.0 ENDIF AA(IDM+1,J) = AA(IDM,J) + A00 ENDDO !j DO I= 1,N IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF AXY = (1.0-DX)*(1.0-DY)*AA(IP, JP) + + (1.0-DX)* DY *AA(IP, JPP1) + + DX *(1.0-DY)*AA(IPP1,JP) + + DX * DY *AA(IPP1,JPP1) PLON(I) = AXY PLON(I) = MOD( PLON(I) + 1080.0, 360.0 ) IF (PLON(I).GT.180.0) THEN PLON(I) = PLON(I) - 360.0 ENDIF IF (IDEBUG.EQ.I) THEN WRITE(6,*) 'debug: i,plon = ',I,PLON(I) ENDIF ENDDO !i C READ(11,REC=2,IOSTAT=IOS) A !plat #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read regional.grid.a' WRITE(6,*) 'ios = ',ios CALL EXIT(5) STOP ENDIF AA(1:IDM,:) = A(:,:) DO J= 1,JDM AA(IDM+1,J) = AA(IDM,J) ENDDO !j DO I= 1,N IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF AXY = (1.0-DX)*(1.0-DY)*AA(IP, JP) + + (1.0-DX)* DY *AA(IP, JPP1) + + DX *(1.0-DY)*AA(IPP1,JP) + + DX * DY *AA(IPP1,JPP1) PLAT(I) = AXY IF (IDEBUG.EQ.I) THEN WRITE(6,*) 'debug: i,plat = ',I,PLAT(I) ENDIF ENDDO !i C READ(11,REC=9,IOSTAT=IOS) A !pang #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read regional.grid.a' WRITE(6,*) 'ios = ',ios CALL EXIT(5) STOP ENDIF DO I= 1,N IP = NINT(XP(I)) JP = NINT(YP(I)) IF (IP.GT.IDM) THEN IP = 1 !periodic wrap ENDIF IF (JP.GT.JDM) THEN JP = JDM ENDIF PANG(I) = A(IP,JP) !no interpolation of angle, use nearest point IF (IDEBUG.EQ.I) THEN WRITE(6,*) 'debug: i,pang = ',I,PANG(I) ENDIF ENDDO !i CLOSE(UNIT=11) DEALLOCATE( AA ) C C OPEN ARCHIVE ".a" FILE. C OPEN(UNIT=11, FILE=CFILEA, FORM='UNFORMATTED', STATUS='OLD', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(CFILEA) WRITE(6,*) 'ios = ',ios WRITE(6,*) 'nrecl = ',nrecl CALL EXIT(3) STOP ENDIF C C SURFACE FIELDS (JUST FSD,UB,VB). C IOS = 0 C IF (NLAY.LE.6) THEN LFATAL = .FALSE. IF (I_ARCH(2).LT.0) THEN READ(11,REC=2,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(2).GT.0) THEN READ(11,REC=I_ARCH(2),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF LANDF = 0 LANDL = 0 DO I= 1,N IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN !1-d interpolation IPP1 = IP ELSEIF (DX.EQ.1.0) THEN !1-d interpolation IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN !1-d interpolation JPP1 = JP ELSEIF (DY.EQ.1.0) THEN !1-d interpolation JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF IF (MAX( A(IP, JP), + A(IP, JPP1), + A(IPP1,JP), + A(IPP1,JPP1) ).GE.SPVAL) THEN C C NEAR LAND, IS NEAREST POINT ACTUALLY LAND? C IF (DX.LE.0.5 .AND. DY.LE.0.5) THEN IP = IP JP = JP ELSEIF (DX.LE.0.5 .AND. DY.GT.0.5) THEN IP = IP JP = JPP1 ELSEIF (DX.GT.0.5 .AND. DY.LE.0.5) THEN IP = IPP1 JP = JP ELSE ! (DX.GT.0.5 .AND. DY.GT.0.5) THEN IP = IPP1 JP = JPP1 ENDIF * IF (.TRUE.) THEN !require 4 good points IF (A(IP,JP).GE.SPVAL) THEN C C NEAREST POINT IS LAND. C FSD(I) = SPVAL IF (LANDF.EQ.0) THEN LANDF = I ENDIF LANDL = I IF (IDEBUG.LT.0) THEN WRITE(6,'(a,2i6)') 'nearest point is land: ',IP,JP LDONE = .FALSE. DO JJ= MAX(1,JP-1),MIN(JDM,JP+1) DO II= MAX(1,IP-1),MIN(IDM,IP+1) IF (A(II,JJ).LT.SPVAL) THEN LDONE = .TRUE. WRITE(6,'(a,2i6)') ' but sea at: ',II,JJ ENDIF ENDDO ENDDO IF (.NOT.LDONE) THEN !25-pt search DO JJ= MAX(1,JP-2),MIN(JDM,JP+2) DO II= MAX(1,IP-2),MIN(IDM,IP+2) IF (A(II,JJ).LT.SPVAL) THEN LDONE = .TRUE. WRITE(6,'(a,2i6)') ' but sea at: ',II,JJ ENDIF ENDDO ENDDO ENDIF ENDIF !IDEBUG ELSE C C USE NEAREST POINT AS THE EXACT VALUE. C IF (LANDF.NE.0) THEN WRITE(6,'(a,i6,a,i6,a,i6,a)') + 'WARNING - STATIONS',LANDF, + ' TO',LANDL, + ' OF',N, ' ARE TREATED AS LAND' LANDF = 0 ENDIF FSD(I) = A(IP,JP)/9.806 XP(I) = IP !for all subsequent interpolations YP(I) = JP !for all subsequent interpolations ENDIF ELSE !all 4 points ok IF (LANDF.NE.0) THEN WRITE(6,'(a,i6,a,i6,a,i6,a)') + 'WARNING - STATIONS',LANDF, + ' TO',LANDL, + ' OF',N, ' ARE TREATED AS LAND' LANDF = 0 ENDIF AXY = (1.0-DX)*(1.0-DY)*A(IP, JP) + + (1.0-DX)* DY *A(IP, JPP1) + + DX *(1.0-DY)*A(IPP1,JP) + + DX * DY *A(IPP1,JPP1) FSD(I) = AXY/9.806 ENDIF IF (IDEBUG.EQ.I) THEN WRITE(6,*) 'debug: i,k,fsd= ',I,K,FSD(I) ENDIF ENDDO !i IF (LANDF.NE.0) THEN WRITE(6,'(a,i6,a,i6,a,i6,a)') + 'WARNING - STATIONS',LANDF, + ' TO',LANDL, + ' OF',N, ' ARE TREATED AS LAND' ENDIF IF (I_ARCH(11).LT.0) THEN READ(11,REC=NSURF-1,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(11).GT.0) THEN READ(11,REC=I_ARCH(11),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (IPP1.EQ.IDM) THEN IPP2 = 1 ELSE IPP2 = IPP1+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF IF (A(IP, JP) .EQ.SPVAL) A(IP, JP) = 0.0 IF (A(IPP1,JP) .EQ.SPVAL) A(IPP1,JP) = 0.0 IF (A(IPP2,JP) .EQ.SPVAL) A(IPP2,JP) = 0.0 IF (A(IP, JPP1).EQ.SPVAL) A(IP, JPP1) = 0.0 IF (A(IPP1,JPP1).EQ.SPVAL) A(IPP1,JPP1) = 0.0 IF (A(IPP2,JPP1).EQ.SPVAL) A(IPP2,JPP1) = 0.0 IF (ARTYPE.GT.0) THEN A00 = 0.5*(A(IP, JP) + A(IP+1,JP) ) A01 = 0.5*(A(IP, JPP1) + A(IP+1,JPP1) ) A10 = 0.5*(A(IPP1,JP) + A(IPP2,JP) ) A11 = 0.5*(A(IPP1,JPP1) + A(IPP2,JPP1) ) ELSE !p-grid A00 = A(IP, JP) A01 = A(IP, JPP1) A10 = A(IPP1,JP) A11 = A(IPP1,JPP1) ENDIF !artype AXY = (1.0-DX)*(1.0-DY)*A00 + + (1.0-DX)* DY *A01 + + DX *(1.0-DY)*A10 + + DX * DY *A11 UB( I) = AXY IF (A00.EQ.0.0 .OR. A01.EQ.0.0 .OR. + A10.EQ.0.0 .OR. A11.EQ.0.0 ) THEN * LFATAL = .TRUE. * WRITE(6,*) 'ERROR - STATION ',I,' HAS BAD U', * + IP,IPP1,IPP2,JP,JPP1 WRITE(6,*) 'WARNING - STATION ',I,' HAS BAD U', + IP,IPP1,IPP2,JP,JPP1 WRITE(6,*) + 'DX,DY = ',DX,DY WRITE(6,*) + 'U.0 = ',A(IP, JP), A(IPP1,JP), A(IPP2,JP) WRITE(6,*) + 'U.1 = ',A(IP, JPP1),A(IPP1,JPP1),A(IPP2,JPP1) ENDIF ENDDO !i C IF (I_ARCH(12).LT.0) THEN READ(11,REC=NSURF, IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(12).GT.0) THEN READ(11,REC=I_ARCH(12),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF IF (JPP1.NE.JDM) THEN JPP2 = JPP1 + 1 ELSE JPP2 = JPP1 ENDIF IF (A(IP, JP) .EQ.SPVAL) A(IP, JP) = 0.0 IF (A(IP, JP+1).EQ.SPVAL) A(IP, JP+1) = 0.0 IF (A(IP, JPP2).EQ.SPVAL) A(IP, JPP2) = 0.0 IF (A(IPP1,JP) .EQ.SPVAL) A(IPP1,JP) = 0.0 IF (A(IPP1,JP+1).EQ.SPVAL) A(IPP1,JP+1) = 0.0 IF (A(IPP1,JPP2).EQ.SPVAL) A(IPP1,JPP2) = 0.0 IF (ARTYPE.GT.0) THEN A00 = 0.5*(A(IP, JP) + A(IP, JP+1) ) A01 = 0.5*(A(IP, JPP1) + A(IP, JPP2) ) A10 = 0.5*(A(IPP1,JP) + A(IPP1,JP+1) ) A11 = 0.5*(A(IPP1,JPP1) + A(IPP1,JPP2) ) ELSE !p-grid A00 = A(IP, JP) A01 = A(IP, JPP1) A10 = A(IPP1,JP) A11 = A(IPP1,JPP1) ENDIF !artype AXY = (1.0-DX)*(1.0-DY)*A00 + + (1.0-DX)* DY *A01 + + DX *(1.0-DY)*A10 + + DX * DY *A11 VB( I) = AXY P(0,I) = 0.0 IF (A00.EQ.0.0 .OR. A01.EQ.0.0 .OR. + A10.EQ.0.0 .OR. A11.EQ.0.0 ) THEN * LFATAL = .TRUE. * WRITE(6,*) 'ERROR - STATION ',I,' HAS BAD V', * + IP,IPP1,JP,JPP1,JPP2 WRITE(6,*) 'WARNING - STATION ',I,' HAS BAD V', + IP,IPP1,JP,JPP1,JPP2 WRITE(6,*) + 'DX,DY = ',DX,DY WRITE(6,*) + 'V.0 = ',A(IP, JP),A(IP, JP+1),A(IP, JPP2) WRITE(6,*) + 'V.1 = ',A(IPP1,JP),A(IPP1,JP+1),A(IPP1,JPP2) ENDIF ENDDO !i IF (LFATAL) THEN CALL EXIT(9) STOP ENDIF ELSE !mean archive DO I= 1,N UB( I) = 0.0 VB( I) = 0.0 P(0,I) = 0.0 ENDDO !i ENDIF C C LAYER FIELDS. C DO K= 1,KDM KREC0 = NSURF+(NLAY+NDIF+NTRC)*(K-1) c c --- u IF (I_ARCH(13).LT.0) THEN READ(11,REC=KREC0+1,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(13).GT.0) THEN READ(11,REC=I_ARCH(13),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (IPP1.EQ.IDM) THEN IPP2 = 1 ELSE IPP2 = IPP1+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF IF (A(IP, JP) .EQ.SPVAL) A(IP, JP) = 0.0 IF (A(IPP1,JP) .EQ.SPVAL) A(IPP1,JP) = 0.0 IF (A(IPP2,JP) .EQ.SPVAL) A(IPP2,JP) = 0.0 IF (A(IP, JPP1).EQ.SPVAL) A(IP, JPP1) = 0.0 IF (A(IPP1,JPP1).EQ.SPVAL) A(IPP1,JPP1) = 0.0 IF (A(IPP2,JPP1).EQ.SPVAL) A(IPP2,JPP1) = 0.0 IF (ARTYPE.GT.0) THEN A00 = 0.5*(A(IP, JP) + A(IP+1, JP) ) A01 = 0.5*(A(IP, JPP1) + A(IP+1, JPP1) ) A10 = 0.5*(A(IPP1,JP) + A(IPP2, JP) ) A11 = 0.5*(A(IPP1,JPP1) + A(IPP2, JPP1) ) ELSE !p-grid A00 = A(IP, JP) A01 = A(IP, JPP1) A10 = A(IPP1,JP) A11 = A(IPP1,JPP1) ENDIF !artype AXY = (1.0-DX)*(1.0-DY)*A00 + + (1.0-DX)* DY *A01 + + DX *(1.0-DY)*A10 + + DX * DY *A11 UK(K,I) = UB(I) + AXY ENDDO !i c c --- v IF (I_ARCH(14).LT.0) THEN READ(11,REC=KREC0+2,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(14).GT.0) THEN READ(11,REC=I_ARCH(14),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF IF (JPP1.NE.JDM) THEN JPP2 = JPP1 + 1 ELSE JPP2 = JPP1 ENDIF IF (A(IP, JP) .EQ.SPVAL) A(IP, JP) = 0.0 IF (A(IP, JP+1).EQ.SPVAL) A(IP, JP+1) = 0.0 IF (A(IP, JPP2).EQ.SPVAL) A(IP, JPP2) = 0.0 IF (A(IPP1,JP) .EQ.SPVAL) A(IPP1,JP) = 0.0 IF (A(IPP1,JP+1).EQ.SPVAL) A(IPP1,JP+1) = 0.0 IF (A(IPP1,JPP2).EQ.SPVAL) A(IPP1,JPP2) = 0.0 IF (ARTYPE.GT.0) THEN A00 = 0.5*(A(IP, JP) + A(IP, JP+1) ) A01 = 0.5*(A(IP, JPP1) + A(IP, JPP2) ) A10 = 0.5*(A(IPP1,JP) + A(IPP1,JP+1) ) A11 = 0.5*(A(IPP1,JPP1) + A(IPP1,JPP2) ) ELSE !p-grid A00 = A(IP, JP) A01 = A(IP, JPP1) A10 = A(IPP1,JP) A11 = A(IPP1,JPP1) ENDIF !artype AXY = (1.0-DX)*(1.0-DY)*A00 + + (1.0-DX)* DY *A01 + + DX *(1.0-DY)*A10 + + DX * DY *A11 VK(K,I) = VB(I) + AXY ENDDO !i c c --- p IF (ABS(ARTYPE).EQ.1) THEN IF (I_ARCH(15).LT.0) THEN READ(11,REC=KREC0+3,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(15).GT.0) THEN READ(11,REC=I_ARCH(15),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF ELSE !mean archive READ(11,REC=KREC0+4,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF AXY = (1.0-DX)*(1.0-DY)*A(IP, JP) + + (1.0-DX)* DY *A(IP, JPP1) + + DX *(1.0-DY)*A(IPP1,JP) + + DX * DY *A(IPP1,JPP1) P(K,I) = P(K-1,I) + AXY/ONEM !interface depth in m IF (IDEBUG.EQ.I) THEN WRITE(6,*) 'debug: i,k,p = ',I,K,P(K,I) ENDIF ENDDO !i c c --- t IF (ABS(ARTYPE).EQ.1) THEN IF (I_ARCH(16).LT.0) THEN READ(11,REC=KREC0+4,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(16).GT.0) THEN READ(11,REC=I_ARCH(16),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF ELSE !mean archive READ(11,REC=KREC0+5,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF AXY = (1.0-DX)*(1.0-DY)*A(IP, JP) + + (1.0-DX)* DY *A(IP, JPP1) + + DX *(1.0-DY)*A(IPP1,JP) + + DX * DY *A(IPP1,JPP1) TK(K,I) = AXY IF (IDEBUG.EQ.I) THEN WRITE(6,*) 'debug: i,k,tk = ',I,K,TK(K,I) ENDIF ENDDO !i c c --- s IF (ABS(ARTYPE).EQ.1) THEN IF (I_ARCH(17).LT.0) THEN READ(11,REC=KREC0+5,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSEIF (I_ARCH(17).GT.0) THEN READ(11,REC=I_ARCH(17),IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ELSE A(:,:) = 0.0 ENDIF ELSE !mean archive READ(11,REC=KREC0+6,IOSTAT=IOS) A #ifdef ENDIAN_IO CALL ENDIAN_SWAP(A,IDM*JDM) #endif ENDIF IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read ',TRIM(CFILEA) CALL EXIT(4) STOP ENDIF DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN CYCLE ENDIF IP = XP(I) JP = YP(I) DX = XP(I) - IP DY = YP(I) - JP IF (DX.EQ.0.0) THEN IPP1 = IP ELSEIF (DX.EQ.1.0) THEN IPP1 = IP+1 IP = IPP1 ELSEIF (IP.EQ.IDM) THEN IPP1 = 1 ELSE IPP1 = IP+1 ENDIF IF (DY.EQ.0.0) THEN JPP1 = JP ELSEIF (DY.EQ.1.0) THEN JPP1 = JP+1 JP = JPP1 ELSE JPP1 = JP+1 ENDIF AXY = (1.0-DX)*(1.0-DY)*A(IP, JP) + + (1.0-DX)* DY *A(IP, JPP1) + + DX *(1.0-DY)*A(IPP1,JP) + + DX * DY *A(IPP1,JPP1) SK(K,I) = AXY ENDDO !i ENDDO !k C C INTERPOLATE ALL PROFILES TO Z C DO I= 1,N IF (FSD(I).EQ.SPVAL) THEN SST( I) = SPVAL SSS( I) = SPVAL DO K= 1,KZ UZ(K,I) = SPVAL VZ(K,I) = SPVAL TZ(K,I) = SPVAL SZ(K,I) = SPVAL ENDDO !k CYCLE ENDIF C IF (BOT.NE.0.0) THEN C C IGNORE LAYERS WITHIN BOT OF THE BOTTOM C DO K= 2,KDM IF (P(K-1,I).GE.P(KDM,I)-BOT) THEN UK(K,I) = UK(K-1,I) VK(K,I) = VK(K-1,I) TK(K,I) = TK(K-1,I) SK(K,I) = SK(K-1,I) ENDIF ENDDO !k ENDIF C SST(I) = TK(1,I) SSS(I) = SK(1,I) C FLAG = SPVAL IF (ITYPE.EQ.-2) THEN CALL LAYER2Z_PPM(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM, & UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG) ELSEIF (ITYPE.EQ.-1) THEN CALL LAYER2Z_PLM(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM, & UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG) ELSEIF (ITYPE.EQ. 0) THEN CALL LAYER2Z_PCM(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM, & UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG) ELSEIF (ITYPE.EQ. 1 .AND. IDEBUG.EQ.I) THEN CALL LAYER2Z_LIN_DEBUG( & UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM, & UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG) ELSEIF (ITYPE.EQ. 1) THEN CALL LAYER2Z_LIN(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM, & UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG) ELSE WRITE(6,"(2a)") + 'Usage: hycom_profile2z_nc', + ' archv.a list.txt z.txt archz.nc [itype [bot]]' WRITE(6,*) 'unsupported itype value' CALL EXIT(1) STOP ENDIF IF (PANG(I).NE.0.0) THEN C C ROTATE TO EASTWARDS,NORTHWARDS C DO K= 1,KZ IF (UZ(K,I).NE.SPVAL) THEN DX = UZ(K,I) DY = VZ(K,I) UZ(K,I) = COS(PANG(I))*DX + SIN(-PANG(I))*DY VZ(K,I) = COS(PANG(I))*DY - SIN(-PANG(I))*DX ENDIF ENDDO !k ENDIF IF (IDEBUG.EQ.I) THEN DO K= 1,KZ WRITE(6,*) 'debug: i,k,tz = ',I,K,TZ(K,I) ENDDO !k ENDIF ENDDO !i C C NETCDF OUTPUT. C ARTYPE = ABS(ARTYPE) CALL PROFOUT(UZ,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' u-veloc.', ! plot name & 'u', ! ncdf name (mersea) & 'eastward_sea_water_velocity', ! ncdf standard_name & 'm/s', ! units & KZ, CFILEP) CALL PROFOUT(VZ,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' v-veloc.', ! plot name & 'v', ! ncdf name (mersea) & 'northward_sea_water_velocity', ! ncdf standard_name & 'm/s', ! units & KZ, CFILEP) IF (KZ.NE.1 .OR. ZZ(1).NE.0.0) THEN CALL PROFOUT(TZ,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' temp ', ! plot name & 'temperature', ! ncdf name & 'sea_water_potential_temperature', ! ncdf standard_name & 'degC', ! units & KZ, CFILEP) CALL PROFOUT(SZ,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' salinity', ! plot name & 'salinity', ! ncdf name & 'sea_water_salinity', ! ncdf standard_name & 'psu', ! units & KZ, CFILEP) ENDIF CALL PROFOUT(FSD,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' SSH ', ! plot name & 'ssh', ! ncdf name (mersea) & 'sea_surface_elevation', ! ncdf standard_name & 'm', ! units & 1, CFILEP) CALL PROFOUT(SST,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' SST ', ! plot name & 'sst', ! ncdf name & 'sea_surface_temperature', ! ncdf standard_name & 'degC', ! units & 1, CFILEP) CALL PROFOUT(SSS,ZZ, PLON,PLAT,DP,N, & ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT, & ' SSS ', ! plot name & 'sss', ! ncdf name & 'sea_surface_salinity', ! ncdf standard_name & 'psu', ! units & 1, CFILEP) END