PROGRAM SWM_TO_GRADS !=== ROGRAM DOCUMENTATION BLOCK ============================= ! . . . . ! PROGRAM: SWM_TO_GRADS INTERPOLATE SWM FIELDS TO LON_LAT GRADS FIELDS ! PRGMMR: ZUPANSKI ORG: CIRA/CSU DATE: 2004-06-25 ! ! ABSTRACT: COMPUTES THE SWM SCALAR VALUE AT LATITUDE-LONGITUDE GRID ! AND WRITES 2-D FIELDS IN GRADS FORMAT ! ! PROGRAM HISTORY LOG: ! 2004-06-25 M.ZUPANSKI - 2-D shallow-water model scalars (no vectors yet) ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE: IBM SP4 !=================================================================== real,parameter :: deg_incr=4.50 !!! average grid-point increment (degrees) integer,parameter :: i_swm=2562 !!! No. of height grid points !---------------------------------------------------------------- !! real,parameter :: deg_incr=2.25 !!! average grid-point increment (degrees) !! integer,parameter :: i_swm=10242 !!! No. of height grid points !---------------------------------------------------------------- integer,parameter :: swm_in=21 !!! input SWM files (1-D) integer,parameter :: ll_out=51 !!! output GRADS files (2-D) integer,parameter :: j_swm=2*i_swm-4 !!! No. of wind grid points integer :: i_lon,j_lat real,parameter :: PI=3.1415927 real :: xlong integer :: i_h !------------------------------------------------------------------- !--- interpolation weights and addresses real,dimension(:,:),allocatable :: w1,w2,w3,w4 integer,dimension(:,:),allocatable :: k1,k2,k3,k4 !-- SW model variables real,dimension(:),allocatable :: & u,v,urot,vrot real,dimension(:),allocatable :: & ulon,ulat real,dimension(:),allocatable :: & hs,hlon,hlat real,dimension(:),allocatable :: & h,psi,chi,relative,div,ke,eta !-- GRADS variables real,dimension(:,:),allocatable :: & h_grads,hs_grads & ,u_grads,v_grads & ,psi_grads,chi_grads & ,relative_grads,div_grads,ke_grads,eta_grads !====================================================== !---------------------------------------------------------- write(*,*) "START FOR SWM MODEL. No OF PTS=",i_swm write(*,*) "deg_incr=",deg_incr !---------------------------------------------------------- i_lon=INT(360./deg_incr + 0.001) j_lat=INT(180./deg_incr + 0.001) + 1 write(*,*) "GRADS (lon-lat) grid dimensions: i,j=",i_lon,j_lat write(*,*) "GRADS (lon-lat) grid increment (degrees) =",deg_incr !--- get input SWM variables (1-D) write(*,*) "SWM: height dim=",i_swm," wind dim=",j_swm allocate(h(1:i_swm)) ! height allocate(psi(1:i_swm)) ! stream function allocate(chi(1:i_swm)) ! velocity potential allocate(relative(1:i_swm)) ! relative vorticity allocate(div(1:i_swm)) ! divergence allocate(ke(1:i_swm)) ! kinetic energy allocate(eta(1:i_swm)) ! absolute vorticity allocate(hs(1:i_swm)) ! sfc height (topography) allocate(hlon(1:i_swm)) ! height longitude (rad) allocate(hlat(1:i_swm)) ! height latitude(rad) allocate(u(1:j_swm)) ! model u-wind allocate(v(1:j_swm)) ! model v-wind allocate(urot(1:j_swm)) ! Earth u-wind (rotated) allocate(vrot(1:j_swm)) ! Earth v-wind (rotated) allocate(ulon(1:j_swm)) ! wind longitude (rad) allocate(ulat(1:j_swm)) ! wind latitude (rad) REWIND swm_in READ(swm_in) h READ(swm_in) u READ(swm_in) v READ(swm_in) psi READ(swm_in) chi READ(swm_in) relative READ(swm_in) div READ(swm_in) ke READ(swm_in) eta READ(swm_in) hs READ(swm_in) hlon READ(swm_in) hlat READ(swm_in) ulon READ(swm_in) ulat READ(swm_in) urot READ(swm_in) vrot write(*,*) "INPUT: h min,max=",minval(h),maxval(h) write(*,*) "INPUT: u min,max=",minval(u),maxval(u) write(*,*) "INPUT: v min,max=",minval(v),maxval(v) write(*,*) "INPUT: psi min,max=",minval(psi),maxval(psi) write(*,*) "INPUT: chi min,max=",minval(chi),maxval(chi) write(*,*) "INPUT: relative min,max=",minval(relative),maxval(relative) write(*,*) "INPUT: div min,max=",minval(div),maxval(div) write(*,*) "INPUT: ke min,max=",minval(ke),maxval(ke) write(*,*) "INPUT: eta min,max=",minval(eta),maxval(eta) write(*,*) "INPUT: hs min,max=",minval(hs),maxval(hs) write(*,*) "INPUT: urot min,max=",minval(urot),maxval(urot) write(*,*) "INPUT: vrot min,max=",minval(vrot),maxval(vrot) write(*,*) "INPUT: hlat min,max=",minval(hlat),maxval(hlat) write(*,*) "INPUT: hlon min,max=",minval(hlon),maxval(hlon) write(*,*) "INPUT: ulat min,max=",minval(ulat),maxval(ulat) write(*,*) "INPUT: ulon min,max=",minval(ulon),maxval(ulon) !-- Calculate (horizontal) interpolation weights allocate(w1(1:i_lon,1:j_lat)) allocate(w2(1:i_lon,1:j_lat)) allocate(w3(1:i_lon,1:j_lat)) allocate(w4(1:i_lon,1:j_lat)) allocate(k1(1:i_lon,1:j_lat)) allocate(k2(1:i_lon,1:j_lat)) allocate(k3(1:i_lon,1:j_lat)) allocate(k4(1:i_lon,1:j_lat)) call weights & (i_lon,j_lat,i_swm,deg_incr,hlon,hlat,w1,w2,w3,w4,k1,k2,k3,k4) !-- Do interpolation from SWM lon-lat grid to GRADS lon-lat grid !!== height points allocate(h_grads(1:i_lon,1:j_lat)) allocate(psi_grads(1:i_lon,1:j_lat)) allocate(chi_grads(1:i_lon,1:j_lat)) allocate(relative_grads(1:i_lon,1:j_lat)) allocate(div_grads(1:i_lon,1:j_lat)) allocate(ke_grads(1:i_lon,1:j_lat)) allocate(eta_grads(1:i_lon,1:j_lat)) allocate(hs_grads(1:i_lon,1:j_lat)) h_grads(:,:)=0.0 psi_grads(:,:)=0.0 chi_grads(:,:)=0.0 relative_grads(:,:)=0.0 div_grads(:,:)=0.0 eta_grads(:,:)=0.0 ke_grads(:,:)=0.0 hs_grads(:,:)=0.0 call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,h,h_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,psi,psi_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,chi,chi_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,relative,relative_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,div,div_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,ke,ke_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,eta,eta_grads) call intpl (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,hs,hs_grads) deallocate(hlon) ; deallocate(hlat) deallocate(h) deallocate(hs) deallocate(psi) deallocate(chi) deallocate(relative) deallocate(div) deallocate(ke) deallocate(eta) write(*,*) "OUTPUT: h min,max=",minval(h_grads),maxval(h_grads) write(*,*) "OUTPUT: hs min,max=",minval(hs_grads),maxval(hs_grads) write(*,*) "OUTPUT: psi min,max=",minval(psi_grads),maxval(psi_grads) write(*,*) "OUTPUT: chi min,max=",minval(chi_grads),maxval(chi_grads) write(*,*) "OUTPUT: relative min,max=",minval(relative_grads),maxval(relative_grads) write(*,*) "OUTPUT: div min,max=",minval(div_grads),maxval(div_grads) write(*,*) "OUTPUT: ke min,max=",minval(ke_grads),maxval(ke_grads) write(*,*) "OUTPUT: eta min,max=",minval(eta_grads),maxval(eta_grads) !!== wind points !-- Calculate (horizontal) interpolation weights w1=0.0 w2=0.0 w3=0.0 w4=0.0 k1=0.0 k2=0.0 k3=0.0 k4=0.0 call weights & (i_lon,j_lat,j_swm,deg_incr,ulon,ulat,w1,w2,w3,w4,k1,k2,k3,k4) deallocate(ulon) ; deallocate(ulat) allocate(u_grads(1:i_lon,1:j_lat)) allocate(v_grads(1:i_lon,1:j_lat)) u_grads(:,:)=0.0 v_grads(:,:)=0.0 call intpl (i_lon,j_lat,j_swm,w1,w2,w3,w4,k1,k2,k3,k4,urot,u_grads) call intpl (i_lon,j_lat,j_swm,w1,w2,w3,w4,k1,k2,k3,k4,vrot,v_grads) deallocate(u) deallocate(v) write(*,*) "OUTPUT: u min,max=",minval(u_grads),maxval(u_grads) write(*,*) "OUTPUT: v min,max=",minval(v_grads),maxval(v_grads) !--- write in GRADS format REWIND ll_out write(ll_out) h_grads write(ll_out) u_grads write(ll_out) v_grads write(ll_out) psi_grads write(ll_out) chi_grads write(ll_out) relative_grads write(ll_out) div_grads write(ll_out) ke_grads write(ll_out) eta_grads write(ll_out) hs_grads deallocate(h_grads) deallocate(u_grads) deallocate(v_grads) deallocate(hs_grads) deallocate(psi_grads) deallocate(chi_grads) deallocate(relative_grads) deallocate(div_grads) deallocate(ke_grads) deallocate(eta_grads) !----------------------------------------------------------------------- STOP END