| 0,0 → 1,1154 |
|---|
| PROGRAM rotate_lalo |
| |
| c ******************************************************************* |
| c * Rotate to a geographical latitude/longitude grid * |
| c * Michael Sprenger / Autumn 2006 * |
| c ******************************************************************* |
| |
| implicit none |
| |
| c ------------------------------------------------------------------- |
| c Declaration of parameters and variables |
| c ------------------------------------------------------------------- |
| |
| c Specification of input parameters |
| real clon,clat,rotation |
| integer nx,ny |
| real dx,dy,xmin,ymin |
| integer nfield |
| character*80 fieldname(100) |
| |
| c Numerical and physical parameters |
| real degrad |
| parameter (degrad=0.0174532925) |
| real omegaerd |
| parameter (omegaerd=7.292e-5 ) |
| real eps |
| parameter (eps=0.001) |
| |
| c Variables for input Z file : height level |
| character*80 in_cfn |
| real in_varmin(4),in_varmax(4),in_stag(4) |
| integer in_vardim(4) |
| real in_mdv |
| integer in_ndim |
| integer in_nx,in_ny,in_nz |
| real in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy |
| integer in_ntimes |
| real in_aklev(500),in_bklev(500) |
| real in_aklay(500),in_bklay(500) |
| real in_time |
| real in_pollon,in_pollat |
| integer in_nvars |
| character*80 in_vnam(100) |
| integer in_idate(5) |
| real,allocatable, dimension (:,:,:) :: in_field3 |
| real,allocatable, dimension (:,:,:) :: in_vect1 |
| real,allocatable, dimension (:,:,:) :: in_vect2 |
| |
| c Variables for output Z file : height level |
| character*80 out_cfn |
| real out_varmin(4),out_varmax(4),out_stag(4) |
| integer out_vardim(4) |
| real out_mdv |
| integer out_ndim |
| integer out_nx,out_ny,out_nz |
| real out_xmin,out_xmax,out_ymin,out_ymax,out_dx,out_dy |
| integer out_ntimes |
| real out_aklev(500),out_bklev(500) |
| real out_aklay(500),out_bklay(500) |
| real out_time |
| real out_pollon,out_pollat |
| integer out_idate(5) |
| real,allocatable, dimension (:,:,:) :: out_field3 |
| real,allocatable, dimension (:,:,:) :: out_vect1 |
| real,allocatable, dimension (:,:,:) :: out_vect2 |
| real,allocatable, dimension (:,:) :: out_lat,out_lon |
| real,allocatable, dimension (:,:) :: out_x,out_y |
| real,allocatable, dimension (:,:) :: out_coriol |
| |
| c Auxiliary variables |
| integer ifield |
| integer i,j,k |
| integer cdfid,cstid |
| character*80 cfn |
| integer stat,ierr,isok |
| real time |
| character*80 varname,cdfname,varname1,varname2 |
| integer idate(5),stdate(5),datar(14) |
| integer rotmode |
| integer i1,i2,i3,i4 |
| integer isvector |
| real lat |
| character*80 name |
| |
| c Externals |
| real sdis |
| external sdis |
| |
| c ------------------------------------------------------------------- |
| c Preparations |
| c ------------------------------------------------------------------- |
| |
| print*,'*********************************************************' |
| print*,'* rotate_lalo *' |
| print*,'*********************************************************' |
| |
| c Read parameter file |
| open(10,file='fort.10') |
| |
| read(10,*) in_cfn |
| read(10,*) out_cfn |
| |
| read(10,*) name,nx |
| if ( trim(name).ne.'GEO_NX' ) stop |
| read(10,*) name,ny |
| if ( trim(name).ne.'GEO_NY' ) stop |
| read(10,*) name,dx |
| if ( trim(name).ne.'GEO_DX' ) stop |
| read(10,*) name,dy |
| if ( trim(name).ne.'GEO_DY' ) stop |
| read(10,*) name,xmin |
| if ( trim(name).ne.'GEO_XMIN') stop |
| read(10,*) name,ymin |
| if ( trim(name).ne.'GEO_YMIN') stop |
| read(10,*) name,clon |
| if ( trim(name).ne.'CLON' ) stop |
| read(10,*) name,clat |
| if ( trim(name).ne.'CLAT' ) stop |
| read(10,*) name,rotation |
| if ( trim(name).ne.'CROT' ) stop |
| |
| read(10,*) nfield |
| do i=1,nfield |
| read(10,*) fieldname(i) |
| enddo |
| close(10) |
| |
| print*,clon,clat,rotation |
| print*,nx,ny,dx,dy,xmin,ymin |
| print*,trim(in_cfn),' -> ',trim(out_cfn) |
| |
| c Get grid description for input Z file : height level |
| call cdfopn(in_cfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,in_nvars,in_vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname=fieldname(1) |
| call check_varok (isok,varname,in_vnam,in_nvars) |
| if (isok.eq.0) goto 998 |
| call getdef(cdfid,varname,in_ndim,in_mdv,in_vardim, |
| > in_varmin,in_varmax,in_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| in_nx =in_vardim(1) |
| in_ny =in_vardim(2) |
| in_xmin=in_varmin(1) |
| in_ymin=in_varmin(2) |
| call getlevs(cstid,in_nz,in_aklev,in_bklev,in_aklay,in_bklay,ierr) |
| call getgrid(cstid,in_dx,in_dy,ierr) |
| in_xmax=in_xmin+real(in_nx-1)*in_dx |
| in_ymax=in_ymin+real(in_ny-1)*in_dy |
| call gettimes(cdfid,in_time,in_ntimes,ierr) |
| call getstart(cstid,in_idate,ierr) |
| call getpole(cstid,in_pollon,in_pollat,ierr) |
| call clscdf(cstid,ierr) |
| call clscdf(cdfid,ierr) |
| |
| c Set grid description for output file : height level |
| out_vardim(1) = nx |
| out_vardim(2) = ny |
| out_vardim(3) = in_nz |
| out_varmin(1) = xmin |
| out_varmin(2) = ymin |
| out_varmin(3) = in_varmin(3) |
| out_varmax(1) = xmin+real(nx-1)*dx |
| out_varmax(2) = ymin+real(ny-1)*dy |
| out_varmax(3) = in_varmax(3) |
| do i=1,in_nz |
| out_aklay(i) = in_aklay(i) |
| out_bklay(i) = in_bklay(i) |
| out_aklev(i) = in_aklev(i) |
| out_bklev(i) = in_bklev(i) |
| enddo |
| out_dx = dx |
| out_dy = dy |
| out_time = in_time |
| out_ntimes = in_ntimes |
| out_ndim = 4 |
| out_mdv = in_mdv |
| out_nx = out_vardim(1) |
| out_ny = out_vardim(2) |
| out_nz = out_vardim(3) |
| out_xmin = out_varmin(1) |
| out_ymin = out_varmin(2) |
| out_pollon = 0. |
| out_pollat = 90. |
| do i=1,5 |
| out_idate(i) = in_idate(i) |
| enddo |
| |
| c Allocate memory for all fields |
| allocate(in_field3(in_nx,in_ny,in_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array in_field3 ***' |
| allocate(in_vect1(in_nx,in_ny,in_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array in_vect1 ***' |
| allocate(in_vect2(in_nx,in_ny,in_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array in_vect2 ***' |
| allocate(out_field3(out_nx,out_ny,out_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_field3 ***' |
| allocate(out_vect1(out_nx,out_ny,out_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_vect1 ***' |
| allocate(out_vect2(out_nx,out_ny,out_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_vect2 ***' |
| allocate(out_lat(out_nx,out_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_lat ***' |
| allocate(out_lon(out_nx,out_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_lon ***' |
| allocate(out_x(out_nx,out_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_x ***' |
| allocate(out_y(out_nx,out_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_y ***' |
| allocate(out_coriol(out_nx,out_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array out_coriol ***' |
| |
| c Create constants file and output file (if requested) |
| datar(1)=out_nx |
| datar(2)=out_ny |
| datar(3)=int(1000.*out_varmax(2)) |
| datar(4)=int(1000.*out_varmin(1)) |
| datar(5)=int(1000.*out_varmin(2)) |
| datar(6)=int(1000.*out_varmax(1)) |
| datar(7)=int(1000.*out_dx) |
| datar(8)=int(1000.*out_dy) |
| datar(9)=out_nz |
| datar(10)=1 |
| datar(11)=1 |
| datar(12)=0 |
| datar(13)=int(1000.*out_pollon) |
| datar(14)=int(1000.*out_pollat) |
| |
| cfn=trim(out_cfn)//'_cst' |
| call wricst(cfn,datar,out_aklev,out_bklev, |
| > out_aklay,out_bklay,out_idate) |
| |
| call crecdf(trim(out_cfn),cdfid,out_varmin,out_varmax, |
| > out_ndim,trim(cfn),ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cdfid,ierr) |
| |
| c ----------------------------------------------------------------- |
| c Loop through all fields - rotate to new grid |
| c ----------------------------------------------------------------- |
| |
| do ifield=1,nfield |
| |
| c Check if scalar or vectorial field (X for scalar, U.V for vector) |
| varname=fieldname(ifield) |
| i1=1 |
| i2=1 |
| i3=0 |
| i4=0 |
| 100 if (varname(i1:i1).eq.' ') then |
| i1=i1+1 |
| goto 100 |
| endif |
| i2=i1 |
| 101 if ((varname(i2:i2).ne.' ').and.(varname(i2:i2).ne.'.')) then |
| i2=i2+1 |
| goto 101 |
| endif |
| if (varname(i2:i2).eq.'.') then |
| i3=i2+1 |
| 102 if (varname(i3:i3).eq.' ') then |
| i3=i3+1 |
| goto 102 |
| endif |
| i4=i3 |
| 104 if (varname(i4:i4).ne.' ') then |
| i4=i4+1 |
| goto 104 |
| endif |
| endif |
| if ((i3.ne.0).and.(i4.ne.0)) then |
| isvector=1 |
| i2=i2-1 |
| varname1=varname(i1:i2) |
| i4=i4-1 |
| varname2=varname(i3:i4) |
| print*,'Rotating vector : ', |
| > trim(varname1),' / ',trim(varname2) |
| else |
| isvector=0 |
| i2=i2-1 |
| varname=varname(i1:i2) |
| print*,'Rotating scalar : ',trim(varname) |
| endif |
| |
| c Rotate a scalar field |
| if (isvector.eq.0) then |
| |
| c Read input field |
| call cdfopn(in_cfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname,in_ndim,in_mdv,in_vardim, |
| > in_varmin,in_varmax,in_stag,ierr) |
| in_nz=in_vardim(3) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,in_time,0,in_field3,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| |
| c Rotate to new coordinate system |
| out_nz=in_nz |
| call getenvir (clon,clat,rotation,0, |
| > in_field3,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv, |
| > out_field3,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,out_nz) |
| |
| c Write output scalar field |
| call cdfwopn(trim(out_cfn),cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| out_vardim(3)=out_nz |
| call putdef(cdfid,varname,4,out_mdv,out_vardim, |
| > out_varmin,out_varmax,out_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| call putdat(cdfid,varname,out_time,0,out_field3,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cdfid,ierr) |
| |
| c Rotate vector field |
| else if (isvector.eq.1) then |
| |
| c Read input vector field |
| call cdfopn(in_cfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname1,in_ndim,in_mdv,in_vardim, |
| > in_varmin,in_varmax,in_stag,ierr) |
| in_nz=in_vardim(3) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname1,in_time,0,in_vect1,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname2,in_ndim,in_mdv,in_vardim, |
| > in_varmin,in_varmax,in_stag,ierr) |
| in_nz=in_vardim(3) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname2,in_time,0,in_vect2,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| |
| c Get new vector component in x direction |
| out_nz=in_nz |
| call getenvir (clon,clat,rotation,1, |
| > in_vect1,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv, |
| > out_field3,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,out_nz) |
| do i=1,out_nx |
| do j=1,out_ny |
| do k=1,out_nz |
| out_vect1(i,j,k)=out_field3(i,j,k) |
| enddo |
| enddo |
| enddo |
| call getenvir (clon,clat,rotation,2, |
| > in_vect2,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv, |
| > out_field3,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,out_nz) |
| do i=1,out_nx |
| do j=1,out_ny |
| do k=1,out_nz |
| if ( (abs(out_vect1 (i,j,k)-out_mdv).gt.eps).and. |
| > (abs(out_field3(i,j,k)-out_mdv).gt.eps) ) then |
| out_vect1(i,j,k)=out_vect1(i,j,k)- |
| > out_field3(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c Get new vector component in y direction |
| out_nz=in_nz |
| call getenvir (clon,clat,rotation,2, |
| > in_vect1,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv, |
| > out_field3,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,out_nz) |
| do i=1,out_nx |
| do j=1,out_ny |
| do k=1,out_nz |
| out_vect2(i,j,k)=out_field3(i,j,k) |
| enddo |
| enddo |
| enddo |
| call getenvir (clon,clat,rotation,1, |
| > in_vect2,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,in_nz,in_mdv, |
| > out_field3,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,out_nz) |
| do i=1,out_nx |
| do j=1,out_ny |
| do k=1,out_nz |
| if ( (abs(out_vect2 (i,j,k)-out_mdv).gt.eps).and. |
| > (abs(out_field3(i,j,k)-out_mdv).gt.eps) ) then |
| out_vect2(i,j,k)=out_vect2(i,j,k)+ |
| > out_field3(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c Write output vector field |
| call cdfwopn(trim(out_cfn),cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| out_vardim(3)=out_nz |
| call putdef(cdfid,varname1,4,out_mdv,out_vardim, |
| > out_varmin,out_varmax,out_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| call putdat(cdfid,varname1,out_time,0,out_vect1,ierr) |
| if (ierr.ne.0) goto 997 |
| call putdef(cdfid,varname2,4,out_mdv,out_vardim, |
| > out_varmin,out_varmax,out_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| call putdat(cdfid,varname2,out_time,0,out_vect2,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cdfid,ierr) |
| |
| endif |
| |
| enddo |
| |
| c ----------------------------------------------------------------- |
| c Write additional fields: latitude, longitude, Coriolis parameter |
| c ----------------------------------------------------------------- |
| |
| c Open the output file |
| call cdfwopn(trim(out_cfn),cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Geographical latitude |
| varname='RLAT' |
| print*,'Rotating scalar : ',trim(varname) |
| do i=1,in_nx |
| do j=1,in_ny |
| in_field3(i,j,1)=in_ymin+real(j-1)*in_dy |
| enddo |
| enddo |
| call getenvir (clon,clat,rotation,0, |
| > in_field3,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,1,in_mdv, |
| > out_lat,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,1) |
| out_vardim(3)=1 |
| call putdef(cdfid,varname,4,out_mdv,out_vardim, |
| > out_varmin,out_varmax,out_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| call putdat(cdfid,varname,out_time,0,out_lat,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Geographical longitude |
| varname='RLON' |
| print*,'Rotating scalar : ',trim(varname) |
| do i=1,in_nx |
| do j=1,in_ny |
| in_field3(i,j,1)=in_xmin+real(i-1)*in_dx |
| enddo |
| enddo |
| call getenvir (clon,clat,rotation,0, |
| > in_field3,in_dx,in_dy, |
| > in_xmin,in_ymin,in_nx,in_ny,1,in_mdv, |
| > out_lon,out_dx,out_dy, |
| > out_xmin,out_ymin,out_nx,out_ny,1) |
| out_vardim(3)=1 |
| call putdef(cdfid,varname,4,out_mdv,out_vardim, |
| > out_varmin,out_varmax,out_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| call putdat(cdfid,varname,out_time,0,out_lon,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| |
| |
| c Close output file |
| call clscdf(cdfid,ierr) |
| |
| c ----------------------------------------------------------------- |
| c Exception handling and format specs |
| c ----------------------------------------------------------------- |
| |
| stop |
| |
| 998 print*,'Z: Problems with input rotated grid' |
| stop |
| |
| 997 print*,'Z: Problems with output to lat/lon grid' |
| stop |
| |
| end |
| |
| |
| c ******************************************************************************** |
| c * SUBROUTINE: ROTATION TO LATITUDE/LONGITUDE COORDINATE SYSTEM * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Subroutine to get environment of strcof |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE getenvir (clon,clat,rotation,rotmode, |
| > inar, rdx,rdy,rxmin,rymin,rnx,rny,rnz,mdv, |
| > outar, dx, dy, xmin, ymin, nx, ny, nz) |
| |
| c Rotate from a local quasi-cartesian coordiante system into lat/lon coordinate |
| c system. |
| |
| implicit none |
| |
| c Declaration of input parameters |
| integer rnx,rny,rnz |
| real rdx,rdy,rxmin,rymin |
| real inar(rnx,rny,rnz) |
| real clon,clat,rotation |
| real mdv |
| integer rotmode |
| |
| c Declaration of output parameters |
| integer nx,ny,nz |
| real dx,dy,xmin,ymin |
| real outar(nx,ny,nz) |
| |
| c Set numerical and physical constants |
| real g2r |
| parameter (g2r=0.0174533) |
| real pi180 |
| parameter (pi180=3.14159265359/180.) |
| real eps |
| parameter (eps=0.0001) |
| |
| |
| c Flag for test mode |
| integer test |
| parameter (test=0) |
| character*80 testfile |
| parameter (testfile='ROTGRID') |
| |
| c Auxiliary variables |
| real pollon,pollat |
| integer i,j,k,l |
| real rlon(nx,ny),rlat(nx,ny) |
| real rlon1(nx,ny),rlat1(nx,ny) |
| real lon(nx,ny),lat(nx,ny) |
| real rotangle1(nx,ny),rotangle2(nx,ny) |
| real rotangle(nx,ny) |
| real sinoutar(nx,ny,nz),cosoutar(nx,ny,nz) |
| real rxmax,rymax |
| real xind,yind,pind |
| real outval |
| integer ix,iy |
| real ax,ay,az,bx,by,bz,zx,zy,zz |
| real clon1,clat1,clon2,clat2 |
| real rindx,rindy |
| integer indx,indy,indr,indu |
| real frac0i,frac0j,frac1i,frac1j |
| character*80 cdfname,varname,cstname |
| real vx1,vx2,vy1,vy2,angle,vx2min |
| integer s |
| |
| c Externals |
| real lmtolms,phtophs |
| external lmtolms,phtophs |
| |
| c Get geographical coordinates |
| do i=1,nx |
| do j=1,ny |
| lon(i,j)=xmin+real(i-1)*dx |
| lat(i,j)=ymin+real(j-1)*dy |
| enddo |
| enddo |
| |
| c First rotation |
| pollon=clon-180. |
| if (pollon.lt.-180.) pollon=pollon+360. |
| pollat=90.-clat |
| do i=1,nx |
| do j=1,ny |
| rlon1(i,j)=lmtolms(lat(i,j),lon(i,j),pollat,pollon) |
| rlat1(i,j)=phtophs(lat(i,j),lon(i,j),pollat,pollon) |
| enddo |
| enddo |
| |
| c Second coordinate transformation |
| pollon=-180. |
| pollat=90.+rotation |
| do i=1,nx |
| do j=1,ny |
| rlon(i,j)=90.+lmtolms(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon) |
| rlat(i,j)=phtophs(rlat1(i,j),rlon1(i,j)-90.,pollat,pollon) |
| enddo |
| enddo |
| |
| c Get the angle between the rotated and non-rotated coordinate frame |
| if ((rotmode.eq.1).or.(rotmode.eq.2)) then |
| do i=2,nx-1 |
| do j=2,ny-1 |
| |
| c Angle between latitude circles |
| vx1=1. |
| vy1=0. |
| vx2min=(rlon(i+1,j)-rlon(i-1,j))*cos(pi180*rlat(i,j)) |
| do s=-1,1,2 |
| vx2=(rlon(i+1,j)-rlon(i-1,j)+real(s)*360.)* |
| > cos(pi180*rlat(i,j)) |
| if (abs(vx2).lt.abs(vx2min)) vx2min=vx2 |
| enddo |
| vx2=vx2min |
| vy2=rlat(i+1,j)-rlat(i-1,j) |
| call getangle(vx1,vy1,vx2,vy2,angle) |
| rotangle1(i,j)=angle |
| |
| c Angle between longitude circles |
| vx1=0. |
| vy1=1. |
| vx2min=(rlon(i,j+1)-rlon(i,j-1))*cos(pi180*rlat(i,j)) |
| do s=-1,1,2 |
| vx2=(rlon(i+1,j)-rlon(i-1,j)+real(s)*360.)* |
| > cos(pi180*rlat(i,j)) |
| if (abs(vx2).lt.abs(vx2min)) vx2min=vx2 |
| enddo |
| vx2=vx2min |
| vy2=rlat(i,j+1)-rlat(i,j-1) |
| call getangle(vx1,vy1,vx2,vy2,angle) |
| rotangle2(i,j)=angle |
| |
| enddo |
| enddo |
| |
| c Set the angle at the boundaries |
| do i=1,nx |
| rotangle1(i,ny)=2.0*rotangle1(i,ny-1)-rotangle1(i,ny-2) |
| rotangle1(i,1) =2.0*rotangle1(i,2)-rotangle1(i,3) |
| rotangle2(i,ny)=2.0*rotangle2(i,ny-1)-rotangle2(i,ny-2) |
| rotangle2(i,1) =2.0*rotangle2(i,2)-rotangle2(i,3) |
| enddo |
| do j=1,ny |
| rotangle1(nx,j)=2.0*rotangle1(nx-1,j)-rotangle1(nx-2,j) |
| rotangle1(1,j) =2.0*rotangle1(2,j)-rotangle1(3,j) |
| rotangle2(nx,j)=2.0*rotangle2(nx-1,j)-rotangle2(nx-2,j) |
| rotangle2(1,j) =2.0*rotangle2(2,j)-rotangle2(3,j) |
| enddo |
| |
| c Set the final rotation angle |
| do i=1,nx |
| do j=1,ny |
| rotangle(i,j)=0.5*(rotangle1(i,j)+rotangle2(i,j)) |
| enddo |
| enddo |
| |
| endif |
| |
| c Bring longitude into domain of geographical grid (shift by 360 deg) |
| do i=1,nx |
| do j=1,ny |
| 100 if (rlon(i,j).lt.rxmin) then |
| rlon(i,j)=rlon(i,j)+360. |
| goto 100 |
| endif |
| 102 if (rlon(i,j).gt.(rxmin+real(rnx-1)*rdx)) then |
| rlon(i,j)=rlon(i,j)-360. |
| goto 102 |
| endif |
| enddo |
| enddo |
| |
| c Rotate the scalar or the vector component |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| |
| c Get index |
| rindx=(rlon(i,j)-rxmin)/rdx+1. |
| rindy=(rlat(i,j)-rymin)/rdy+1. |
| indx=int(rindx) |
| indy=int(rindy) |
| if ((indx.lt.1).or.(indx.gt.rnx).or. |
| > (indy.lt.1).or.(indy.gt.rny)) then |
| outar(i,j,k)=mdv |
| goto 103 |
| endif |
| |
| c Get inidices of left and upper neighbours |
| indr=indx+1 |
| if (indr.gt.rnx) indr=1 |
| indu=indy+1 |
| if (indu.gt.rny) indu=ny |
| |
| c Do linear interpolation |
| if ( ( abs(inar(indx ,indy, k)-mdv).gt.eps ).and. |
| & ( abs(inar(indx ,indu, k)-mdv).gt.eps ).and. |
| & ( abs(inar(indr ,indy, k)-mdv).gt.eps ).and. |
| & ( abs(inar(indr ,indu, k)-mdv).gt.eps ) ) then |
| frac0i=rindx-float(indx) |
| frac0j=rindy-float(indy) |
| frac1i=1.-frac0i |
| frac1j=1.-frac0j |
| outval = inar(indx ,indy, k ) * frac1i * frac1j |
| & + inar(indx ,indu, k ) * frac1i * frac0j |
| & + inar(indr ,indy, k ) * frac0i * frac1j |
| & + inar(indr ,indu, k ) * frac0i * frac0j |
| else |
| outval=mdv |
| endif |
| |
| c Update output array |
| outar(i,j,k)=outval |
| |
| c Next |
| 103 continue |
| |
| enddo |
| enddo |
| enddo |
| |
| c Get components for tests |
| if (test.eq.1) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| cosoutar(i,j,k)=outar(i,j,k)*cos(pi180*rotangle(i,j)) |
| sinoutar(i,j,k)=-outar(i,j,k)*sin(pi180*rotangle(i,j)) |
| enddo |
| enddo |
| enddo |
| endif |
| |
| c Get correct component of rotated field |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ( abs(outar(i,j,k)-mdv).gt.eps ) then |
| if (rotmode.eq.1) then |
| outar(i,j,k)= outar(i,j,k)*cos(pi180*rotangle(i,j)) |
| else if (rotmode.eq.2) then |
| outar(i,j,k)=-outar(i,j,k)*sin(pi180*rotangle(i,j)) |
| else if (rotmode.eq.0) then |
| outar(i,j,k)=outar(i,j,k) |
| endif |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c Test mode: Write grids to cdf file |
| if (test.eq.1) then |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='RLON1' |
| call writecdf2D(cdfname,cstname, |
| > varname,rlon1,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 1,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='RLON' |
| call writecdf2D(cdfname,cstname, |
| > varname,rlon,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='LON' |
| call writecdf2D(cdfname,cstname, |
| > varname,lon,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='RLAT1' |
| call writecdf2D(cdfname,cstname, |
| > varname,rlat1,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='RLAT' |
| call writecdf2D(cdfname,cstname, |
| > varname,rlat,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='LAT' |
| call writecdf2D(cdfname,cstname, |
| > varname,lat,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='ANGLE1' |
| call writecdf2D(cdfname,cstname, |
| > varname,rotangle1,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='ANGLE2' |
| call writecdf2D(cdfname,cstname, |
| > varname,rotangle2,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='U' |
| call writecdf2D(cdfname,cstname, |
| > varname,cosoutar,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| cdfname=testfile |
| cstname=trim(testfile)//'_cst' |
| varname='V' |
| call writecdf2D(cdfname,cstname, |
| > varname,sinoutar,0., |
| > rdx,rdy,rxmin,rymin,rnx,rny, |
| > 0,1) |
| |
| endif |
| |
| END |
| |
| |
| c -------------------------------------------------------------------------------- |
| c Auxiliary routines: angle between two vectors |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE getangle (vx1,vy1,vx2,vy2,angle) |
| |
| c Given two vectors <vx1,vy1> and <vx2,vy2>, determine the angle (in deg) |
| c between the two vectors. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| real vx1,vy1 |
| real vx2,vy2 |
| real angle |
| |
| c Auxiliary variables and parameters |
| real len1,len2,len3 |
| real val1,val2,val3 |
| real pi |
| parameter (pi=3.14159265359) |
| |
| len1=sqrt(vx1*vx1+vy1*vy1) |
| len2=sqrt(vx2*vx2+vy2*vy2) |
| |
| if ((len1.gt.0.).and.(len2.gt.0.)) then |
| vx1=vx1/len1 |
| vy1=vy1/len1 |
| vx2=vx2/len2 |
| vy2=vy2/len2 |
| |
| val1=vx1*vx2+vy1*vy2 |
| val2=-vy1*vx2+vx1*vy2 |
| |
| len3=sqrt(val1*val1+val2*val2) |
| |
| if ( (val1.ge.0.).and.(val2.ge.0.) ) then |
| val3=acos(val1/len3) |
| else if ( (val1.lt.0.).and.(val2.ge.0.) ) then |
| val3=pi-acos(abs(val1)/len3) |
| else if ( (val1.ge.0.).and.(val2.le.0.) ) then |
| val3=-acos(val1/len3) |
| else if ( (val1.lt.0.).and.(val2.le.0.) ) then |
| val3=-pi+acos(abs(val1)/len3) |
| endif |
| else |
| val3=0. |
| endif |
| |
| angle=180./pi*val3 |
| |
| END |
| |
| c -------------------------------------------------------------------------------- |
| c Transformation routine: LMSTOLM and PHSTOPH from library gm2em |
| c -------------------------------------------------------------------------------- |
| |
| REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM) |
| C |
| C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H% |
| C |
| C**** LMTOLMS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM |
| C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) |
| C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT |
| C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM) |
| C** AUFRUF : LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM) |
| C** ENTRIES : KEINE |
| C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF |
| C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM |
| C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT |
| C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM) |
| C** VERSIONS- |
| C** DATUM : 03.05.90 |
| C** |
| C** EXTERNALS: KEINE |
| C** EINGABE- |
| C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM |
| C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM |
| C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS |
| C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS |
| C** AUSGABE- |
| C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION |
| C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0) |
| C** |
| C** COMMON- |
| C** BLOECKE : KEINE |
| C** |
| C** FEHLERBE- |
| C** HANDLUNG : KEINE |
| C** VERFASSER: G. DE MORSIER |
| |
| REAL LAM,PHI,POLPHI,POLLAM |
| |
| DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 / |
| |
| ZSINPOL = SIN(ZPIR18*POLPHI) |
| ZCOSPOL = COS(ZPIR18*POLPHI) |
| ZLAMPOL = ZPIR18*POLLAM |
| ZPHI = ZPIR18*PHI |
| ZLAM = LAM |
| IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0 |
| ZLAM = ZPIR18*ZLAM |
| |
| ZARG1 = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI) |
| ZARG2 = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI) |
| IF (ABS(ZARG2).LT.1.E-30) THEN |
| IF (ABS(ZARG1).LT.1.E-30) THEN |
| LMTOLMS = 0.0 |
| ELSEIF (ZARG1.GT.0.) THEN |
| LMTOLMS = 90.0 |
| ELSE |
| LMTOLMS = -90.0 |
| ENDIF |
| ELSE |
| LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2) |
| ENDIF |
| |
| RETURN |
| END |
| |
| |
| REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM) |
| C |
| C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H% |
| C |
| C**** PHTOPHS - FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI |
| C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) |
| C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT |
| C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM) |
| C** AUFRUF : PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM) |
| C** ENTRIES : KEINE |
| C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF |
| C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM |
| C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT |
| C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM) |
| C** VERSIONS- |
| C** DATUM : 03.05.90 |
| C** |
| C** EXTERNALS: KEINE |
| C** EINGABE- |
| C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM |
| C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM |
| C** POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS |
| C** POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS |
| C** AUSGABE- |
| C** PARAMETER: ROTIERTE BREITE PHIS ALS WERT DER FUNKTION |
| C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0) |
| C** |
| C** COMMON- |
| C** BLOECKE : KEINE |
| C** |
| C** FEHLERBE- |
| C** HANDLUNG : KEINE |
| C** VERFASSER: G. DE MORSIER |
| |
| REAL LAM,PHI,POLPHI,POLLAM |
| |
| DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 / |
| |
| ZSINPOL = SIN(ZPIR18*POLPHI) |
| ZCOSPOL = COS(ZPIR18*POLPHI) |
| ZLAMPOL = ZPIR18*POLLAM |
| ZPHI = ZPIR18*PHI |
| ZLAM = LAM |
| IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0 |
| ZLAM = ZPIR18*ZLAM |
| ZARG = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI) |
| |
| PHTOPHS = ZRPI18*ASIN(ZARG) |
| |
| RETURN |
| END |
| |
| |
| c --------------------------------------------------------- |
| c Spherical distance between two lat/lon points |
| c --------------------------------------------------------- |
| |
| real function sdis(xp,yp,xq,yq) |
| c |
| c calculates spherical distance (in km) between two points given |
| c by their spherical coordinates (xp,yp) and (xq,yq), respectively. |
| c |
| real re |
| parameter (re=6370.) |
| real xp,yp,xq,yq,arg |
| |
| arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq) |
| if (arg.lt.-1.) arg=-1. |
| if (arg.gt.1.) arg=1. |
| sdis=re*acos(arg) |
| |
| end |
| |
| |
| c ******************************************************************************** |
| c * NETCDF INPUT/OUTPUT * |
| c ******************************************************************************** |
| |
| c --------------------------------------------------------- |
| c Subroutines to write the netcdf output file |
| c --------------------------------------------------------- |
| |
| subroutine writecdf2D(cdfname,cstname, |
| > varname,arr,ctime, |
| > dx,dy,xmin,ymin,nx,ny, |
| > crefile,crevar) |
| |
| c Create and write to the netcdf file <cdfname>. The variable |
| c with name <varname> and with time <time> is written. The data |
| c are in the two-dimensional array <arr>. The list <dx,dy,xmin, |
| c ymin,nx,ny> specifies the output grid. The flags <crefile> and |
| c <crevar> determine whether the file and/or the variable should |
| c be created. 1: create / 0: not created |
| |
| IMPLICIT NONE |
| |
| c Declaration of input parameters |
| character*80 cdfname,cstname,varname |
| integer nx,ny |
| real arr(nx,ny) |
| real dx,dy,xmin,ymin |
| integer ctime |
| integer crefile,crevar |
| |
| c Further variables |
| real varmin(4),varmax(4),stag(4) |
| integer ierr,cdfid,ndim,vardim(4) |
| real mdv |
| integer datar(14),stdate(5) |
| integer i |
| real time |
| |
| C Definitions |
| varmin(1)=xmin |
| varmin(2)=ymin |
| varmin(3)=1050. |
| varmax(1)=xmin+real(nx-1)*dx |
| varmax(2)=ymin+real(ny-1)*dy |
| varmax(3)=1050. |
| ndim=4 |
| vardim(1)=nx |
| vardim(2)=ny |
| vardim(3)=1 |
| stag(1)=0. |
| stag(2)=0. |
| stag(3)=0. |
| mdv=-999.98999 |
| time=real(ctime) |
| |
| C Create the file |
| if (crefile.eq.0) then |
| call cdfwopn(cdfname,cdfid,ierr) |
| if (ierr.ne.0) goto 906 |
| else if (crefile.eq.1) then |
| call crecdf(cdfname,cdfid, |
| > varmin,varmax,ndim,cstname,ierr) |
| if (ierr.ne.0) goto 903 |
| |
| C Write the constants file |
| datar(1)=vardim(1) |
| datar(2)=vardim(2) |
| datar(3)=int(1000.*varmax(2)) |
| datar(4)=int(1000.*varmin(1)) |
| datar(5)=int(1000.*varmin(2)) |
| datar(6)=int(1000.*varmax(1)) |
| datar(7)=int(1000.*dx) |
| datar(8)=int(1000.*dy) |
| datar(9)=1 |
| datar(10)=1 |
| datar(11)=0 ! data version |
| datar(12)=0 ! cstfile version |
| datar(13)=0 ! longitude of pole |
| datar(14)=90000 ! latitude of pole |
| do i=1,5 |
| stdate(i)=0 |
| enddo |
| call wricst(cstname,datar,0.,0.,0.,0.,stdate) |
| endif |
| |
| c Write the data to the netcdf file, and close the file |
| if (crevar.eq.1) then |
| call putdef(cdfid,varname,ndim,mdv, |
| > vardim,varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 904 |
| endif |
| call putdat(cdfid,varname,time,0,arr,ierr) |
| if (ierr.ne.0) goto 905 |
| call clscdf(cdfid,ierr) |
| |
| return |
| |
| c Error handling |
| 903 print*,'*** Problem to create netcdf file ***' |
| stop |
| 904 print*,'*** Problem to write definition ***' |
| stop |
| 905 print*,'*** Problem to write data ***' |
| stop |
| 906 print*,'*** Problem to open netcdf file ***' |
| stop |
| |
| END |
| |
| c ---------------------------------------------------------------- |
| c Check whether variable is found on netcdf file |
| c ---------------------------------------------------------------- |
| |
| subroutine check_varok (isok,varname,varlist,nvars) |
| |
| c Check whether the variable <varname> is in the list <varlist(nvars)>. |
| c If this is the case, <isok> is incremented by 1. Otherwise <isok> |
| c keeps its value. |
| |
| implicit none |
| |
| c Declaraion of subroutine parameters |
| integer isok |
| integer nvars |
| character*80 varname |
| character*80 varlist(nvars) |
| |
| c Auxiliary variables |
| integer i |
| |
| c Main |
| do i=1,nvars |
| if (trim(varname).eq.trim(varlist(i))) isok=isok+1 |
| enddo |
| |
| end |
| Property changes: |
|---|
| Added: svn:executable |
|---|