MODULE datashare ! Main variables INTEGER, PARAMETER:: nx = 301, ny = 301, ns = 29, nh = 41 ! Input variables from TCLAPS outmdl*.nc REAL, DIMENSION (nx,ny,ns) :: us, vs, ws, ts, hs, qvs, qis, qws, qrs, qss, qgs REAL :: sig(ns) ! Common variables for both files REAL, DIMENSION (nx) :: rlon REAL, DIMENSION (ny) :: rlat REAL, DIMENSION (nx,ny):: ps, sst,hlowest ! Output variables for H*.nc ! ws is omega (pa/s) on sigma level ! wh is omega (pa/s) on h level ! wwh is W (m/s) on h level ! hlowest is the height of the model lowest sigma level, just (in case) for reference REAL, DIMENSION (nx,ny,nh) :: uh, vh, wh, th, qvh, qih, qwh, qrh, qsh, qgh, ph,wwh REAL, DIMENSION (nx,ny,nh) :: rvor, avor,pvor,thae REAL :: hh(nh) !******************************** REAL, DIMENSION (nx,ny,ns) :: psig,tv INTEGER ncid, IDvar END MODULE datashare ! !--------------------------------------------------------------------------- ! MODULE axi INTEGER, PARAMETER:: nr = 81 INTEGER:: IDnx, IDny, IDnz,IDnt INTEGER:: IDxvar, IDyvar, IDzvar,IDtvar INTEGER:: IDvtm, IDvtsd,IDvrm,IDvrsd,IDwm,IDwsd,IDrvm,IDrvsd,IDavm,IDavsd,IDpvm,IDpvsd ! Azimuthal averages REAL, DIMENSION (nr,nh) :: Axmean, Axsd, Axtemp INTEGER, DIMENSION (nx,ny) :: mask,dist,var2d INTEGER, DIMENSION (nr) :: npoint REAL, DIMENSIOn (nr):: rr END MODULE axi ! ! ************************************************************************** PROGRAM sig2ht ! ! -------------------------------------------------------------------------- ! ! PURPOSE: ! Convert TCLAPS model output from sigma to z coordinate ! Method: linear interpolation in height ! Input fields: ! SST, ps, U, V, W, T, Qv, Qi, Qw, Qr, Qs, Qg ! Output fields: ! SST, ps, Ph, Uh, Vh, Wh, Th, Qvh, Qih, Qrh, Qsh, Qgh ! Mai started to code on 13/12/2008, based on net2ctl.f90 ! ! -------------------------------------------------------------------------- USE datashare USE axi include 'netcdf.inc' ! Input and output file names CHARACTER(50) :: flnin, flnout,flntrk,flnrzout INTEGER :: valdate, valtime REAL :: RdonG, undef,rtemp,dr,dh RdonG = 287.05/9.8 undef = 1.e32 dh = 250. dr = 5.e3 ! *********************************************************************************** ! Get input parmaters ! ------------------- flnin = 'fsig.nc' flnout = 'fhght.nc' flntrk = 'track.trk' flnrzout = 'axirz.nc' Call getvars(flnin) ! ************************************************************************************ ! CONVERT from sigma to height coordinates !----------------- CALL s2h(dh) ! OUTPUT TO netcdf file of the interpolated 3D variables CALL w2Dh(flnout) ! Calculate diagnostic variables CALL calRVor() CALL calAvor() CALL calPVor() CALL calThae() ! Azimuthal averages CALL azim(dr) ! Setup 2D netcdf files Call wrzinit(fnlrzout,IDfout) Call axistat(vth,axmean, axsd) ret = NF_PUT_VAR_REAL(IDfout,IDvtm,Axmean(1,1)) ret = NF_PUT_VAR_REAL(IDfout,IDvtsd,Axsd(1,1)) Call axistat(vrh,axmean, axsd) ret = NF_PUT_VAR_REAL(IDfout,IDvrm,Axmean(1,1)) ret = NF_PUT_VAR_REAL(IDfout,IDvrsd,Axsd(1,1)) Call axistat(wh,axmean, axsd) ret = NF_PUT_VAR_REAL(IDfout,IDwm,Axmean(1,1)) ret = NF_PUT_VAR_REAL(IDfout,IDwsd,Axsd(1,1)) Call axistat(wh,axmean, axsd) ret = NF_PUT_VAR_REAL(IDfout,IDwm,Axmean(1,1)) ret = NF_PUT_VAR_REAL(IDfout,IDwsd,Axsd(1,1)) STOP END SUBROUTINE calThae() USE datashare, only: nx,ny,nh,th,qvh,ph,thae REAL:: gg, cappa, hloncp REAL, DIMENSION (nx,ny,nh):: thv gg = 9.8 cappa=0.288 hloncp = 2501./1.0046 thv(:,:,:) = ( th(:,:,:)*(1.+qvh(:,:,:)/.622))*(1.e5/ph(:,:,:))**cappa thae(:,:,:) = thv(:,:,:)*exp(hloncp*qvh(:,:,:)/th(:,:,:)) RETURN END SUBROUTINE calThae SUBROUTINE wrzinit(flnrzout,IDfout) USE axi USE datashare, only:: nh,hh CHARACTER(50)::flnrzout include 'netcdf.inc' ! Create output file iret = NF_CREATE (flnrzout,nf_clobber,IDfout) if (iret.ne.0) then write(*,*) 'Error in opening file ',flnrzout,' iret = ',iret endif ! Put into define mode iret=NF_REDEF(IDfout) write(*,*) 'going to define mode, NF_REDEF, iret = ', iret ! Define dimensions iret = NF_DEF_DIM(IDfout,'x',nr,IDnx) ! iret = NF_DEF_DIM(IDfout,'y',1,IDny) iret = NF_DEF_DIM(IDfout,'z',nh,IDnz) iret = NF_DEF_DIM(IDfout,'t',1,IDnt) ! Define coordinate variables iret = NF_DEF_VAR(IDfout,'x',nf_float,1,IDnx,IDxvar) ! iret = NF_DEF_VAR(IDfout,'y',nf_float,0,1,IDyvar) iret = NF_DEF_VAR(IDfout,'z',nf_float,1,IDnz,IDzvar) iret = NF_DEF_VAR(IDfout,'t',nf_int,0,1,IDtvar) ! Define additional variables iret=NF_DEF_VAR(IDfout,'vtm',nf_float,3,(/IDnx,IDnz,IDnt/),IDvtm) iret=NF_DEF_VAR(IDfout,'vtsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDvtsd) iret=NF_DEF_VAR(IDfout,'vrm',nf_float,3,(/IDnx,IDnz,IDnt/),IDvrm) iret=NF_DEF_VAR(IDfout,'vrsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDvrsd) iret=NF_DEF_VAR(IDfout,'wm',nf_float,3,(/IDnx,IDnz,IDnt/),IDwm) iret=NF_DEF_VAR(IDfout,'wsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDwsd) iret=NF_DEF_VAR(IDfout,'rvm',nf_float,3,(/IDnx,IDnz,IDnt/),IDrvm) iret=NF_DEF_VAR(IDfout,'rvsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDrvsd) iret=NF_DEF_VAR(IDfout,'avm',nf_float,3,(/IDnx,IDnz,IDnt/),IDavm) iret=NF_DEF_VAR(IDfout,'avsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDavsd) iret=NF_DEF_VAR(IDfout,'pvm',nf_float,3,(/IDnx,IDnz,IDnt/),IDpvm) iret=NF_DEF_VAR(IDfout,'pvsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDpvsd) iret=NF_DEF_VAR(IDfout,'vtm',nf_float,3,(/IDnx,IDnz,IDnt/),IDvtm) iret=NF_DEF_VAR(IDfout,'vtsd',nf_float,3,(/IDnx,IDnz,IDnt/),IDvtsd) ! Get out of the define mode iret=NF_ENDDEF(IDfout) ! Put dimension variables in iret = NF_PUT_VAR_REAL(IDfout,IDxvar,rr(1)) iret = NF_PUT_VAR_REAL(IDfout,IDzvar,hh(1)) RETURN END SUBROUTINE wrzinit SUBROUTINE w2Dh(flnout) USE datashare include 'netcdf.inc' CHARACTER(50):: flnout INTEGER:: IDnx, IDny, IDnz,IDnt, IDxvar,IDyvar,IDzvar, IDtvar, & IDp,IDu,IDv,IDomega,IDw,IDt,IDqv,IDqi,IDqw,IDqr,IDqs,IDqg, & INIDsst,IDps,IDhlo ! Create output file iret = NF_CREATE (flnout,nf_clobber,IDfout) if (iret.ne.0) then write(*,*) 'Error in opening file ',flnout,' iret = ',iret endif ! Put into define mode iret=NF_REDEF(IDfout) write(*,*) 'going to define mode, NF_REDEF, iret = ', iret ! Define dimensions iret = NF_DEF_DIM(IDfout,'x',nx,IDnx) iret = NF_DEF_DIM(IDfout,'y',ny,IDny) iret = NF_DEF_DIM(IDfout,'z',nh,IDnz) iret = NF_DEF_DIM(IDfout,'t',1,IDnt) ! Define coordinate variables iret = NF_DEF_VAR(IDfout,'x',nf_float,1,IDnx,IDxvar) iret = NF_DEF_VAR(IDfout,'y',nf_float,1,IDny,IDyvar) iret = NF_DEF_VAR(IDfout,'z',nf_float,1,IDnz,IDzvar) iret = NF_DEF_VAR(IDfout,'t',nf_int,0,1,IDtvar) ! Define additional variables iret=NF_DEF_VAR(IDfout,'ph',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDp) iret=NF_DEF_VAR(IDfout,'uh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDu) iret=NF_DEF_VAR(IDfout,'vh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDv) iret=NF_DEF_VAR(IDfout,'omegah',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDomega) iret=NF_DEF_VAR(IDfout,'wh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDw) iret=NF_DEF_VAR(IDfout,'th',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDt) iret=NF_DEF_VAR(IDfout,'qvh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDqv) iret=NF_DEF_VAR(IDfout,'qih',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDqi) iret=NF_DEF_VAR(IDfout,'qwh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDqw) iret=NF_DEF_VAR(IDfout,'qrh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDqr) iret=NF_DEF_VAR(IDfout,'qsh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDqs) iret=NF_DEF_VAR(IDfout,'qgh',nf_float,4,(/IDnx,IDny,IDnz,IDnt/),IDqg) iret=NF_DEF_VAR(IDfout,'sst',nf_float,3,(/IDnx,IDny,IDnt/),IDsst) iret=NF_DEF_VAR(IDfout,'ps',nf_float,3,(/IDnx,IDny,IDnt/),IDps) iret=NF_DEF_VAR(IDfout,'hlo',nf_float,3,(/IDnx,IDny,IDnt/),IDhlo) ! Get out of DEFINE mode iret=NF_ENDDEF(IDfout) write(*,*) 'Leaving define mode, NF_ENDEF, iret = ', iret ! Put values out iret = NF_PUT_VAR_REAL(IDfout,IDxvar,rlon(1)) iret = NF_PUT_VAR_REAL(IDfout,IDyvar,rlat(1)) iret = NF_PUT_VAR_REAL(IDfout,IDzvar,hh(1)) iret = NF_PUT_VAR_REAL(IDfout,IDp,ph(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDu,uh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDv,vh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDomega,wh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDw,wwh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDt,th(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDqv,qvh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDqi,qih(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDqw,qwh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDqr,qrh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDqs,qsh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDqg,qgh(1,1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDsst,sst(1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDps,ps(1,1)) iret = NF_PUT_VAR_REAL(IDfout,IDhlo,hlowest(1,1)) ! Close output file iret = NF_CLOSE(IDfout) write(*,*) 'Done with file ',flnout RETURN END SUBROUTINE w2Dh SUBROUTINE s2h(dh) USE datashare !include 'netcdf.inc' ! REAL :: RdonG, undef,rtemp RdonG = 287.05/9.8 undef = 1.e32 ! Assign vertical coordinates hh do ih = 1,nh hh(ih) = (ih-1)*dh enddo ! Calculate pressure and Tv at sigma levels do is = 1, ns psig(:,:,is) = sig(is)*ps(:,:) tv(:,:,is) = ts(:,:,is)*(1. + 0.61*qvs(:,:,is)) enddo ! 1D interpolation do i=1,nx do j=1,ny ! Calculate height at sigma levels ! For the bottom model level (ns) dlnp=log(ps(i,j))-log(psig(i,j,ns)) hs(i,j,ns)=RdonG*dlnp*(tv(i,j,ns)+sst(i,j))/2. do ks = ns -1,1,-1 kp = ks +1 dlnp=log(psig(i,j,kp))-log(psig(i,j,ks)) hs(i,j,ks) = hs(i,j,kp) + 0.5*(RdonG*dlnp*(tv(i,j,kp)+tv(i,j,ks))) enddo ! Vertical loop in height do kh = 1,nh ! Find layer where a particular h level stays levlo = 0 do ks = ns,1,-1 if (hh(kh).eq.hs(i,j,ks)) then levlo=ks whi = 0. exit endif ksup = ks -1 rtemp = (hh(kh)-hs(i,j,ks))*(hh(kh)-hs(i,j,ksup)) if (rtemp.lt.0.) then levlo = ks whi = (hh(kh)-hs(i,j,ks))/(hs(i,j,ksup)-hs(i,j,ks)) exit endif enddo ! Lelel found, do interpolation wlo = 1.-whi if (levlo.ne.0) then levhi = levlo -1 uh(i,j,kh) = wlo*us(i,j,levlo) + whi*us(i,j,levhi) vh(i,j,kh) = wlo*vs(i,j,levlo) + whi*vs(i,j,levhi) wh(i,j,kh) = wlo*ws(i,j,levlo) + whi*ws(i,j,levhi) ph(i,j,kh) = wlo*psig(i,j,levlo) + whi*psig(i,j,levhi) th(i,j,kh) = wlo*ts(i,j,levlo) + whi*ts(i,j,levhi) qvh(i,j,kh) = wlo*qvs(i,j,levlo) + whi*qvs(i,j,levhi) qih(i,j,kh) = wlo*qis(i,j,levlo) + whi*qis(i,j,levhi) qwh(i,j,kh) = wlo*qws(i,j,levlo) + whi*qws(i,j,levhi) qrh(i,j,kh) = wlo*qrs(i,j,levlo) + whi*qrs(i,j,levhi) qsh(i,j,kh) = wlo*qss(i,j,levlo) + whi*qss(i,j,levhi) qgh(i,j,kh) = wlo*qgs(i,j,levlo) + whi*qgs(i,j,levhi) ! Convert omega (pa/s) to W(m/s) ! W=-(R/g)*Tv/P * omega ! wwh(i,j,kh) = -RdonG*th(i,j,kh)*(1.+0.61*qvh(i,j,kh))/ph(i,j,kh)*wh(i,j,kh) else uh(i,j,kh) = undef vh(i,j,kh) = undef wh(i,j,kh) = undef wwh(i,j,kh) = undef ph(i,j,kh) = undef th(i,j,kh) = undef qvh(i,j,kh) = undef qih(i,j,khrea) = undef qwh(i,j,kh) = undef qrh(i,j,kh) = undef qsh(i,j,kh) = undef qgh(i,j,kh) = undef endif enddo ! Completed interpolation for each grid point (1D array) uh(i,j,1) = us(i,j,ns) vh(i,j,1) = vs(i,j,ns) wh(i,j,1) = ws(i,j,ns) wwh(i,j,1) = 0. ph(i,j,1) = ps(i,j) hlowest(i,j) = hs(i,j,ns) th(i,j,1) = ts(i,j,ns) qvh(i,j,1) = qvs(i,j,ns) qih(i,j,1) = qis(i,j,ns) qwh(i,j,1) = qws(i,j,ns) qrh(i,j,1) = qrs(i,j,ns) qsh(i,j,1) = qss(i,j,ns) qgh(i,j,1) = qgs(i,j,ns) enddo enddo ! *************************************************************************** RETURN END SUBROUTINE s2h SUBROUTINE getvars() USE datashare include 'netcdf.inc' iret = NF_OPEN (flnin,nf_read, ncid) IF ( iret .ne. 0 ) then write(*,*) '**** ERROR in netCDF2ctl : IRET = ', iret write(*,*) ' Opening NC file ', flnin, ' using NF_OPEN' ENDIF WRITE (*,*) 'Open netCDF file ', ncid ! Inquire information iret= NF_INQ (ncid, ndims, nvars, ngatts, unlimid) IF ( iret .ne. 0 ) then write(*,*) '**** ERROR in netCDF2ctl : IRET = ', iret index = -1 ENDIF ! Get xdimension iret=NF_INQ_VARID (ncid,'lon',IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, rlon(1)) ! Get ydimension iret=NF_INQ_VARID (ncid,'lat',IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, rlat(1)) ! Get zdimension iret=NF_INQ_VARID (ncid,'lvl',IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, sig(1)) ! Get valid_date iret=NF_INQ_VARID (ncid,'valide_date', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, valdate) ! Get valid_time iret=NF_INQ_VARID (ncid,'valide_time', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, valtime) ! Get U iret=NF_INQ_VARID (ncid,'zonal_wnd', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, us(1,1,1)) ! Get V iret=NF_INQ_VARID (ncid,'merid_wnd', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, vs(1,1,1)) ! Get omega iret=NF_INQ_VARID (ncid,'omega', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, ws(1,1,1)) ! Get temperature iret=NF_INQ_VARID (ncid,'air_temp', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, ts(1,1,1)) ! Get vapour mixing ratio iret=NF_INQ_VARID (ncid,'mix_rto', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, qvs(1,1,1)) ! Get qi iret=NF_INQ_VARID (ncid,'cld_ice', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, qis(1,1,1)) ! Get qw iret=NF_INQ_VARID (ncid,'cld_water', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, qws(1,1,1)) ! Get qr iret=NF_INQ_VARID (ncid,'cld_rain', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, qrs(1,1,1)) ! Get qs iret=NF_INQ_VARID (ncid,'cld_snow', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, qss(1,1,1)) ! Get qg iret=NF_INQ_VARID (ncid,'cld_grpl', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, qgs(1,1,1)) ! Get ps iret=NF_INQ_VARID (ncid,'sfc_pres', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, ps(1,1)) ! Get sst iret=NF_INQ_VARID (ncid,'sfc_temp', IDvar) iret=NF_GET_VAR_REAL(ncid, IDvar, sst(1,1)) ! Close input file iret = NF_CLOSE(ncid) RETURN END SUBROUTINE getvars SUBROUTINE azim(dr,flntrk) USE datashare, only:: rr,dist,mask,npoint,nx,ny,nr CHARACTER(50):: flntrk REAL:: Re, fac Re = 6377.79e3 fac = 3.14/180. ! Setting up the azimuthal grid Do i=1,nr rr(i)= (i-1)*dr Enddo ! Find center open(unit= 30,file=flntrk,access='READ',status = 'OLD') Do read (30,'(I3,x,I8,I3,x,A12,2F6.1,F7.1)') nn,idate,ihour,ttx,rlat,rlon,rpmin If (idate == valdate and ihour == valtime) write (*,*) valdate, valtime,rlat,rlon,rpmin close (30) EXIT endif Enddo tclat = rlat tclon = rlon ! Calculate distance do i=1,nx do j=1,ny rtemp = ((rlon(i)-tclon)*cos(fac*rlat(j)))**2 + (rlat(j)-tclat)**2 dist(i,h) = Re*fac*sqrt(rtemp) enddo enddo ! Mapping npoint = 0 do ir = 1,nr ri = rr(ir) ro = rr(ir+1) where (dist >= ri .AND. dist < ro) mask = ir enddo do ir=1,nr do i = 1, nx do j = 1, ny if (mask(i,j).eq.ir) then npoint(ir) = npoint(ir) +1 endif enddo endo enddo RETURN END SUBROUTINE azim