module met use const use comm use var implicit real*8 (a-h,o-z) # include "pnetcdf.inc" contains subroutine metgen dimension fg(j2t),tanphig(j2t), & yg(j0t),yvg(j0t),yvdegg(j0t),ydegg(j0t), & csg(j0t),csvg(j1t),ocsg(j0t),ocsvg(j1t), & dxg(j0t),dxvg(j1t),odxg(j0t),odxvg(j1t), & dyg(j0t),dyvg(j1t),odyg(j0t),odyvg(j1t) real*8 dd,ddv,flag !! kilmer ! tlz= depth of bottom z-level ! ztop= depth of interface between top two layers ! dydx= ratio of dy to dx at each latitude (constant) ! dxmnut= longitude resolution (minutes) ! r0= earth's radius (cm) ! dxdeg= longitudinal resolution (degrees) data tlz_depth,d,ztop,dydx,dxmnut,y0deg,r0/ & 5800.E2,.04, 20.E2, 1., 60.,-86.7,6.4E8/ ! kilmer dxmnut=360.*60./dble(i2t) dxdeg=dxmnut/60. !dydx=dydx/dxdeg westdeg=360./dble(ngx0) westdeg=westdeg*dble(mylon) if (westdeg .ge. 360.) westdeg=westdeg-360. ! open(9,file='data/westdeg'//grd) ! rewind 9 ! write(9,'(f9.4)') westdeg ! close(9) call wr8nc_westdeg !dog ! open(18,file='data/latitudes'//grd,form='unformatted') ! open(10,file='data/yzgrid'//grd,form='unformatted') ! rewind 10 ! rewind 18 tmp=.01*tlz_depth dzz=tmp/(2.*dble(k1)) coeff_c=.001 coeff_d=d coeff_a=tmp*(1.-coeff_d)/(1.-exp(coeff_c*tmp)) coeff_b=-coeff_a zz=0. do k=1,k0+k1 z(k)=100.*(coeff_a+coeff_b*exp(coeff_c*zz)+coeff_d*zz) zz=zz+dzz enddo if(myid.eq.0)write(*,'(10f8.0)') z ! calculate vertical metrics do k=1,k1 odz(k)=1./(z(2*k+1)-z(2*k-1)) enddo ! odzw(1,i,j) and odzw(k0,i,j) are not used do k=2,k1 odzw(k)=1./(z(2*k)-z(2*k-2)) enddo it=mylon*i2 jt=mylat*j2 ! degrees to radians pi_180=3.141592654/180. yvdegg(1)=y0deg ydegg(1)=y0deg yvg(1)=0. dyg(1)=0. odyg(1)=0. ! equatorial dx dx0=r0*dxdeg*pi_180 ! kilmer !dd=21.0/100.0 !ddv=dd/100.0 !flag=1.0 ! kilmer end do j=2,j0t ! initial guess ! yvdegg(j)=yvdegg(j-1)+dxdeg*dydx*cos(yvdegg(j-1)*pi_180) yvdegg(j)=yvdegg(j-1)+1.476-dxdeg*dydx*cos(yvdegg(j-1)*pi_180) ! kilmer ! print*,j,yvdegg(j) dxg(j)=dx0*cos(.5*(yvdegg(j-1)+yvdegg(j))*pi_180) tmp=yvdegg(j) ! yvdegg(j)=yvdegg(j-1)+dydx*dxg(j)/(r0*pi_180) yvdegg(j)=yvdegg(j-1)+1.476-dydx*dxg(j)/(r0*pi_180) ! kilmer ydegg(j)=.5*(yvdegg(j-1)+yvdegg(j)) csg(j)=cos(ydegg(j)*pi_180) dxrad=dxdeg*pi_180*csg(j) dxg(j)=dxrad*r0 dyrad=(yvdegg(j)-yvdegg(j-1))*pi_180 temp=dyrad/dxrad dydeg=yvdegg(j)-yvdegg(j-1) dyg(j)=dydeg*r0*pi_180 dyvg(j-1)=(ydegg(j)-ydegg(j-1))*r0*pi_180 odyg(j)=1./dyg(j) odyvg(j-1)=1./dyvg(j-1) yvg(j)=yvg(j-1)+dyg(j) yg(j)=.5*(yvg(j)+yvg(j-1)) ocsg(j)=1./csg(j) odxg(j)=1./dxg(j) csvg(j-1)=cos(pi_180*yvdegg(j-1)) ocsvg(j-1)=1./csvg(j-1) dxvg(j-1)=dx0*csvg(j-1) odxvg(j-1)=1./dxvg(j-1) ! kilmer !odxvg(j-1)=1./dxvg(j-1) !if (j>j0t-100) then ! dd=dd+ddv ! flag=0.0 !endif !if (flag==1.0) dd=dd-ddv !if (dd<0.0) dd=0.0 enddo ! kilmer end y1deg=yvdegg(j1) jt=mylat*j2 do j=1,j1 csv(j)=csvg(jt+j) ocsv(j)=ocsvg(jt+j) dxv(j)=dxvg(jt+j) odxv(j)=odxvg(jt+j) dyv(j)=dyvg(jt+j) odyv(j)=odyvg(jt+j) enddo do j=1,j0 y(j)=yg(jt+j) yv(j)=yvg(jt+j) yvdeg(j)=yvdegg(jt+j) ydeg(j)=ydegg(jt+j) cs(j)=csg(jt+j) ocs(j)=ocsg(jt+j) dx(j)=dxg(jt+j) odx(j)=odxg(jt+j) dy(j)=dyg(jt+j) ody(j)=odyg(jt+j) enddo ! open(16,file='data/cell_centers'//grd) if (myid ==0) write(*,16) (j,ydeg(j),j=2,j1) 16 format(' j ','cell center latitude'/(i3,f11.4)) ! write(18) ydeg ! ---------------------------- ! set coriolis parameter array ! ---------------------------- omega2=3.141592654/21600. do j=2,j1 ! spherical curvature parameter tanphi(j-1)=tan(pi_180*ydeg(j))/r0 f(j-1)=omega2*sin(pi_180*ydeg(j)) enddo do i=1,i0 xdeg(i)=westdeg+(i-1.5)*(dxmnut/60.) !w xdeg(i)=(i-1.5)*(dxmnut/60.) if(xdeg(i).ge.360.)then xdeg(i)=xdeg(i)-360. elseif(xdeg(i).lt.0.)then xdeg(i)=xdeg(i)+360. endif enddo write(*,'(10f9.4)') (xdeg(i),i=2,i1) ! write(10) j0,k0 write(*,*) j0,k0 ! write(10) tlz_depth,z,odz,odzw,dxmnut,xdeg,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv call wr8nc_yzgrid(tlz_depth,dxmnut) !dog ! open(11,file='data/ydeg'//grd,form='unformatted') ! write(11) ydeg call wr8nc_ydeg !dog ! kilmer !if (myid==0) then ! open(3333,file='data/global_ydegg.log',form='formatted') ! write(3333,322) (ydegg(j),j=1,j0t) ! close(3333) !endif !322 format(/'global ydegg'/(10f10.4)) ! kilmer end if(mylon.eq.0)then write(*,20) (yvdeg(j),j=1,j1) 20 format(/'model latitude grid lines'/(10f8.4)) write(*,21) f 21 format('coriolis parameter'/(1p,10e9.2)) write(*,99) (j,yvdeg(j),j=1,j1) 99 format('j,yvdeg(j)=',i3,f9.4) endif end subroutine subroutine wr8nc_westdeg use const, only:ngx0,ngy0 use comm, only: mylon,myid,m_cart implicit none integer :: rank,ncid,ierr,nout INTEGER(KIND=MPI_OFFSET_KIND) :: start(1),count(1),bufcount=1,dimm integer :: x_dim integer :: dims=1,dimids(1),varid real*8 :: mydata(1) dimm = ngx0*ngy0 mydata(1)=360./real(ngx0) mydata(1)=mydata(1)*real(mylon) if (mydata(1) > 360.) mydata(1)=mydata(1)-360. ierr=nfmpi_create(m_cart,'data/westdeg.nc',nf_64bit_offset,mpi_info_null,ncid) nout = nfmpi_def_dim(ncid,"myid+1",dimm,x_dim) dimids = (/x_dim/) nout = nfmpi_def_var(ncid,"westdeg",NF_DOUBLE,dims,dimids,varid) nout = nfmpi_enddef(ncid) start = (/myid+1/) count = (/1/) nout = nfmpi_put_vara_all(ncid,varid,start,count,mydata,bufcount,MPI_DOUBLE) nout = nfmpi_close(ncid) end subroutine wr8nc_westdeg subroutine wr8nc_ydeg use const, only: ngy0,j0 use comm, only: myid,m_cart implicit none character :: fn*12 integer :: ncid,ierr,ns integer :: j0id integer :: iydeg integer(kind=mpi_offset_kind) :: dim_nn integer(kind=mpi_offset_kind) :: ns_i,nc_i !output file name fn = 'data/ydeg.nc' !create ncfile and id ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid) !define dimension and assign dimID dim_nn = j0*ngy0 ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn,j0id) !define variable ns = j0id ierr=nfmpi_def_var(ncid,"ydeg",nf_double,1,ns,iydeg) ierr=nfmpi_enddef(ncid) !put variable ns_i = mylat*j0+1 nc_i = j0 ierr=nfmpi_put_vara_double_all(ncid,iydeg,ns_i,nc_i,ydeg) ierr=nfmpi_close(ncid) end subroutine wr8nc_ydeg subroutine wr8nc_yzgrid(tlz_depth,dxmnut) use const, only: i0,j0,j1,k0,k1,ngx0,ngy0 use comm, only: m_cart,mylat,mylon implicit none !declare pnetcdf var character :: fn*14 integer :: ierr,ncid,ns integer :: siid,k0id,k1id,kkid,i0id,j0id,j1id integer :: itlz,idxm,iz,iodz,iodzw,ixdeg,iij,iik integer :: iy,iyv,iyvd,iydeg,ics,iocs,icsv,iocsv,idx,iodx integer :: idxv,iodxv,idy,idyv,iody,iodyv integer(kind=mpi_offset_kind) :: dim_nn(7) integer(kind=mpi_offset_kind) :: ns_i,nc_i real*8 :: tlz_depth,dxmnut integer :: tmp(1) !output file name fn = 'data/yzgrid.nc' !create ncfile and id ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid) !define dimension and assign dimID dim_nn(1) = 1 dim_nn(2) = k0 dim_nn(3) = k1 dim_nn(4) = k0+k1 dim_nn(5) = i0*ngx0 dim_nn(6) = j0*ngy0 dim_nn(7) = j1*ngy0 ierr=nfmpi_def_dim(ncid,"pt_val",dim_nn(1),siid) ierr=nfmpi_def_dim(ncid,"k0",dim_nn(2),k0id) ierr=nfmpi_def_dim(ncid,"k1",dim_nn(3),k1id) ierr=nfmpi_def_dim(ncid,"k0+k1",dim_nn(4),kkid) ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(5),i0id) ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(6),j0id) ierr=nfmpi_def_dim(ncid,"j1_global",dim_nn(7),j1id) !define variables and assign varID ns = siid ierr=nfmpi_def_var(ncid,"j0",nf_int,1,ns,iij) ierr=nfmpi_def_var(ncid,"k0",nf_int,1,ns,iik) ierr=nfmpi_def_var(ncid,"tlz_depth",nf_double,1,ns,itlz) ierr=nfmpi_def_var(ncid,"dxmnut",nf_double,1,ns,idxm) ns = kkid ierr=nfmpi_def_var(ncid,"z",nf_double,1,ns,iz) ns = k1id ierr=nfmpi_def_var(ncid,"odz",nf_double,1,ns,iodz) ns = k0id ierr=nfmpi_def_var(ncid,"odzw",nf_double,1,ns,iodzw) ns = i0id ierr=nfmpi_def_var(ncid,"xdeg",nf_double,1,ns,ixdeg) ns = j0id ierr=nfmpi_def_var(ncid,"y",nf_double,1,ns,iy) ierr=nfmpi_def_var(ncid,"yv",nf_double,1,ns,iyv) ierr=nfmpi_def_var(ncid,"yvdeg",nf_double,1,ns,iyvd) ierr=nfmpi_def_var(ncid,"ydeg",nf_double,1,ns,iydeg) ierr=nfmpi_def_var(ncid,"cs",nf_double,1,ns,ics) ierr=nfmpi_def_var(ncid,"ocs",nf_double,1,ns,iocs) ierr=nfmpi_def_var(ncid,"dx",nf_double,1,ns,idx) ierr=nfmpi_def_var(ncid,"odx",nf_double,1,ns,iodx) ierr=nfmpi_def_var(ncid,"dy",nf_double,1,ns,idy) ierr=nfmpi_def_var(ncid,"ody",nf_double,1,ns,iody) ns = j1id ierr=nfmpi_def_var(ncid,"csv",nf_double,1,ns,icsv) ierr=nfmpi_def_var(ncid,"ocsv",nf_double,1,ns,iocsv) ierr=nfmpi_def_var(ncid,"dxv",nf_double,1,ns,idxv) ierr=nfmpi_def_var(ncid,"odxv",nf_double,1,ns,iodxv) ierr=nfmpi_def_var(ncid,"dyv",nf_double,1,ns,idyv) ierr=nfmpi_def_var(ncid,"odyv",nf_double,1,ns,iodyv) ierr=nfmpi_enddef(ncid) !put variables ns_i = 1 nc_i = 1 ierr=nfmpi_begin_indep_data(ncid) tmp(1) = j0 ierr=nfmpi_put_var_int(ncid,iij,tmp) tmp(1) = k0 ierr=nfmpi_put_var_int(ncid,iik,tmp) ierr=nfmpi_put_var_double(ncid,itlz,tlz_depth) ierr=nfmpi_put_var_double(ncid,idxm,dxmnut) nc_i = k0+k1 ierr=nfmpi_put_var_double(ncid,iz,z) nc_i = k1 ierr=nfmpi_put_var_double(ncid,iodz,odz) nc_i = k0 ierr=nfmpi_put_var_double(ncid,iodzw,odzw) ierr=nfmpi_end_indep_data(ncid) ns_i = mylon*i0+1 nc_i = i0 ierr=nfmpi_put_vara_double_all(ncid,ixdeg,ns_i,nc_i,xdeg) ns_i = mylat*j0+1 nc_i = j0 ierr=nfmpi_put_vara_double_all(ncid,iy,ns_i,nc_i,y) ierr=nfmpi_put_vara_double_all(ncid,iyv,ns_i,nc_i,yv) ierr=nfmpi_put_vara_double_all(ncid,iyvd,ns_i,nc_i,yvdeg) ierr=nfmpi_put_vara_double_all(ncid,iydeg,ns_i,nc_i,ydeg) ierr=nfmpi_put_vara_double_all(ncid,ics,ns_i,nc_i,cs) ierr=nfmpi_put_vara_double_all(ncid,iocs,ns_i,nc_i,ocs) ierr=nfmpi_put_vara_double_all(ncid,idx,ns_i,nc_i,dx) ierr=nfmpi_put_vara_double_all(ncid,iodx,ns_i,nc_i,odx) ierr=nfmpi_put_vara_double_all(ncid,idy,ns_i,nc_i,dy) ierr=nfmpi_put_vara_double_all(ncid,iody,ns_i,nc_i,ody) ns_i = mylat*j1+1 nc_i = j1 ierr=nfmpi_put_vara_double_all(ncid,icsv,ns_i,nc_i,csv) ierr=nfmpi_put_vara_double_all(ncid,iocsv,ns_i,nc_i,ocsv) ierr=nfmpi_put_vara_double_all(ncid,idxv,ns_i,nc_i,dxv) ierr=nfmpi_put_vara_double_all(ncid,iodxv,ns_i,nc_i,odxv) ierr=nfmpi_put_vara_double_all(ncid,idyv,ns_i,nc_i,dyv) ierr=nfmpi_put_vara_double_all(ncid,iodyv,ns_i,nc_i,odyv) ierr=nfmpi_close(ncid) end subroutine wr8nc_yzgrid end module