;************************************************* ; roms_5.ncl ; ; Concepts illustrated: ; - Plotting ROMS data ; - Loading NCL functions from another script ; - Changing the reference annotation string for vectors ; - Overlaying vectors and filled contours on a map ; ; roms depth slice using roms_3d_interp with velocity overlay ; note: roms_3d_interp is working on "rho" coords ; if you specify variable "u" or "v" it is transfering them ; first to the "rho" and than interpolate. ; ;************************************************* load "/home/arnold/software/ncl_ncarg/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "/home/arnold/software/ncl_ncarg/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "/home/arnold/software/ncl_ncarg/lib/ncarg/nclscripts/ROMS_utils.ncl" load "/home/arnold/software/ncl_ncarg/lib/ncarg/nclscripts/roms_depth_interp_ratio.ncl" begin ;************************************************ ; User settings ;************************************************ fhis = "/home/arnold/ROMS_AGRIF_MLCS/Roms_tools/Run/FORECAST/roms_his_forecast.nc" ;*********************************************** system("/bin/rm -f /home/arnold/ROMS_AGRIF_MLCS/Roms_tools/Run/FORECAST/roms_interpolation2.nc") ; remove any pre-existing file ncdf = addfile("/home/arnold/ROMS_AGRIF_MLCS/Roms_tools/Run/FORECAST/roms_interpolation2.nc" ,"c") ; open output netCDF file ;=================================================================== ; create global attributes of the file (optional) ;=================================================================== fAtt = True ; assign file attributes fAtt@title = "NCL Simple Approach to netCDF Creation" fAtt@source_file = "roms_his_forecast.nc" fAtt@Conventions = "None" fAtt@creation_date = systemfunc ("date") fileattdef( ncdf, fAtt ) ; copy file attributes ;*********************************************** ; Read data ;*********************************************** his = addfile (fhis,"r") lon2d = his->lon_rho lat2d = his->lat_rho times = his->time ZETA = his->zeta h = his->h hc = his->hc Cs_r = his->Cs_r Sc_r = his->s_rho s = his->salt if(isfilevar(his,"Vtransform")) vtransform = his->Vtransform else print("Do not have Vtransform in file needed for depth calculations, using 2 as default") vtransform = 2 return end if lev = his->s_rho lev_v = (/1.0, 0.75, 0.5, 0.25, 0.0/) T = his->temp salt = his->salt dim_T=dimsizes(T) U1 = his->u U=u2rho(U1) V1 = his->v V = v2rho(V1) ntim = dimsizes(times) ; get dimension sizes nlev = dimsizes(lev_v) nlat1 = dimsizes(lat2d) nlat = nlat1(0) nlon1 = dimsizes(lon2d) nlon = nlon1(1) dimSizes = (/ ntim, nlev, nlat, nlon /) dimSizes2 = (/ ntim, nlat, nlon /) U_out = new(dimSizes,typeof(U)) V_out = new(dimSizes,typeof(V)) T_out = new(dimSizes,typeof(T)) IST_out = new(dimSizes,typeof(T)) S_out = new(dimSizes,typeof(T)) P_out = new(dimSizes,typeof(T)) ;---------- depth calculation ---------------------------------- depth = new(dim_T,typeof(U)) hinv=1./h N=dimSizes(1) M=dimSizes(0) if (vtransform.eq.2) do it=0,M-1 do k=0,N-1 cff = 1./(hc + h); cffr = hc*Sc_r(k) + h*Cs_r(k); depth(it,k,:,:)=ZETA(it,:,:) + ( ZETA(it,:,:) + h )*cffr*cff end do end do end if if (vtransform.eq.1) do it=0,M-1 do k=0,N-1 cffr = hc*(Sc_r(k) - Cs_r(k)) depth(it,k,:,:)=cffr+Cs_r(k)*h + ZETA(it,:,:)*(1+(cffr+Cs_r(k)*h)*hinv) end do end do end if ;-------------- Depth calculation ------------------------------ ;-----------in situ temp calculation --------------------------- ;pre = depth_to_pres(depth,False)*10 pre=abs(depth)*0.01*10 pre@units = "dbar" pref = 0.0 ; user specified reference pressure (dbars) pref@units = "dbar" ; must be decibar units ;istemp = new(dim_T,typeof(T)) ;istemp = potmp_insitu_ocn(T,s,pre, pref, 1 ,True) ; pot(ntim,klev,nlat,nlon) ;------------- in situ temp calculation --------------------------- T_out(:,nlev-1,:,:) = T(:,dim_T(1)-1,:,:) U_out(:,nlev-1,:,:) = U(:,dim_T(1)-1,:,:) V_out(:,nlev-1,:,:) = V(:,dim_T(1)-1,:,:) S_out(:,nlev-1,:,:) = salt(:,dim_T(1)-1,:,:) IST_out(:,nlev-1,:,:) = T(:,dim_T(1)-1,:,:) IST_out@units = "Degrees C" IST_out@standard_name = "sea_water_potential_temperature" IST_out@long_name = "In Situ Temperature" ;IST_out@_FillValue = 9.96921e+36 do it = 0,ntim-1 ; TIME LOOP ; print("Working on time: " + times(it) ) do il = 0,nlev-2 do i=0,dimSizes(2)-1 do j=0,dimSizes(3)-1 ; out = roms_depth_interp_ratio(his,"temp",it,lev_v(il)) ; printVarSummary(out) x=depth(it,:,i,j) min_x = min(x)*lev_v(il) T_out(it,il,i,j) = linint1(x,T(it,:,i,j),False,min_x,0); S_out(it,il,i,j) = linint1(x,salt(it,:,i,j),False,min_x,0); U_out(it,il,i,j) = linint1(x,U(it,:,i,j),False,min_x,0); V_out(it,il,i,j) = linint1(x,V(it,:,i,j),False,min_x,0); P_out(it,il,i,j) = linint1(x,pre(it,:,i,j),False,min_x,0); IST_out(it,il,i,j)= potmp_insitu_ocn(T_out(it,il,i,j),S_out(it,il,i,j),P_out(it,il,i,j), pref, 0 ,True) ; T_out(it,il,:,:) = (/roms_depth_interp_ratio(his,"temp",it,lev_v(il))/) ; S_out(it,il,:,:) = (/roms_depth_interp_ratio(his,"salt",it,lev_v(il))/) ; U_out(it,il,:,:) = (/roms_depth_interp_ratio(his, "u", it,lev_v(il))/) ; they are automatically put to "rho" ; V_out(it,il,:,:) = (/roms_depth_interp_ratio(his, "v", it,lev_v(il))/) end do end do end do end do U_out!0="time" U_out!1="lev" U_out!2="lat" U_out!3="lon" V_out!0="time" V_out!1="lev" V_out!2="lat" V_out!3="lon" T_out!0="time" T_out!1="lev" T_out!2="lat" T_out!3="lon" IST_out!0="time" IST_out!1="lev" IST_out!2="lat" IST_out!3="lon" ; printVarSummary(T_out) ; printVarSummary(U_out) ;===================================================================== ncdf->T = T_out ncdf->U = U_out ncdf->V = V_out ncdf->Time = times ncdf->lat_rho = lat2d ncdf->lon_rho = lon2d ncdf->lev = lev_v ncdf->zeta = ZETA ncdf->istemp = IST_out end