| ! Copyright by http://scriptdemo.blogspot.com MODULE mod_longcheck USE mod_precision IMPLICIT NONE CONTAINS FUNCTION Check180(mylon) ! check whether it is close to the International Date Line (IDL) logical :: Check180 real(kind=rprec) :: mylon if ( abs(abs(mylon)-180.0) < 15.0) then ! here the threshold value is arbitrary Check180=.TRUE. else Check180=.FALSE. endif END FUNCTION Check180 FUNCTION CheckZero(mylon) ! check if whether it is close to the Primary Line logical :: CheckZero real(kind=rprec) :: mylon if ( abs(mylon) < 5.0) then ! here the threshold value is arbitrary CheckZero=.TRUE. else if (abs(abs(mylon)-360.0)< 5.) then CheckZero=.TRUE. else CheckZero=.FALSE. endif endif END FUNCTION CheckZero FUNCTION Back180(mylon) ! convert a longitude into a range of [-180 180] real(kind=rprec) :: mylon real(kind=rprec) :: Back180 Back180=mylon if (Back180<0.0) then do while (Back180 < -180.0) Back180=Back180+360.0 enddo else do while (Back180>180.0) Back180=Back180-360.0 enddo endif return END FUNCTION Back180 FUNCTION Back360(mylon) ! convert a longitude into a range of [0 360] real(kind=rprec) :: mylon real(kind=rprec) :: Back360 Back360=mylon if (Back360<0.0) then do while (Back360 < 0.0) Back360=Back360+360.0 enddo else do while (Back360>360.0) Back360=Back360-360.0 enddo endif return END FUNCTION Back360 SUBROUTINE ll2xy(mylon,mylat,myx,myy) ! convert given long/latitude into x/y in stereographic projection real(kind=rprec) :: mylon,mylat real(kind=rprec) :: myx,myy ! local variables real(kind=rprec) :: mypi,deg2r,myphi,mylam real(kind=rprec) :: rho,Re,phi0,lam0 real(kind=rprec) :: a,b,k,f,e,e2,e4,e6,t Re=6371000.0 ! the earth radius mypi=acos(-1.0) deg2r=mypi/180.0 myphi=mylat*deg2r mylam=mylon*deg2r phi0=0. lam0=0. if (mylat>84.0) then rho=2.0*Re*(sqrt((1.0-sin(myphi))/(1.0+sin(myphi)))) myx=rho*sin(mylam-lam0) myy=-rho*cos(mylam-lam0) myx=myx+2000000.0 myy=myy+2000000.0 elseif (mylat<-84.0) then ! south-hemisphere stereographic projection myphi=-myphi mylam=-mylam lam0=-lam0 rho=2.0*Re*(sqrt((1.0-sin(myphi))/(1.0+sin(myphi)))) myx=rho*sin(mylam-lam0) myy=-rho*cos(mylam-lam0) myx=-myx; myy=-myy myx=myx+2000000.0 myy=myy+2000000.0 else stop "NOT DEFINED YET, maybe just not necessary" endif END SUBROUTINE ll2xy SUBROUTINE xy2ll(myx,myy,mylon,mylat) ! convert given x/y in the stereographic projection into long/latitude real(kind=rprec) :: myx,myy real(kind=rprec) :: mylon,mylat ! local variables real(kind=rprec) :: mypi,deg2r,myphi,mylam,lat0 real(kind=rprec) :: xx,yy real(kind=rprec) :: rho,Re,phi0,lam0 real(kind=rprec) :: a,b,k,f,e,e2,e4,e6,t mypi=acos(1.0) deg2r=mypi/180.0 Re=6371000.0 lam0=0.0 phi0=0.0 lat0=90.0 if (lat0>84.0) then xx=myx-2000000.0 yy=myy-2000000.0 rho = sqrt (xx*xx + yy*yy) myphi=mypi*0.5-2*atan(0.5*rho/Re) mylam=lam0+atan2(xx,-yy) elseif (lat0<-84.0) then stop 'Should not see this msg' xx=-(myx-2000000.0) yy=-(myy-2000000.0) rho=sqrt(xx*xx+yy*yy) myphi=0.5*mypi-2*atan(0.5*rho/Re) myphi=-myphi mylam=-(-lam0+atan2(xx,-yy)) else stop 'Should not see this msg' endif END SUBROUTINE xy2ll END MODULE mod_longcheck |
Thursday, September 30, 2010
Fortran Code dealing with longitude
Saturday, September 25, 2010
Pack netcdf file
| #!/bin/bash # A simple usage of ncpdq from nco (http://nco.sourceforge.net) # Please check more details from the nco website. # Usage: packnc.sh myncfile(s) # NOTE: There will be a loss of precision in type conversion # (converted to short by default) # Copyright by http://scriptdemo.blogspot.com/ IsNCPDQ=`which ncpdq` if [ ${#IsNCPDQ} -eq 0 ]; then echo "ncpdq can not be found in the path..." echo "Please double check your NetCDF Operator (NCO) Path" exit fi if [ $# -eq 0 ]; then sed -n '4p' packnc.sh exit else Nfiles=$* for ncfile in ${Nfiles} do [ -f ${ncfile} ] && ncpdq -P all_xst ${ncfile} tmpout_${ncfile} && mv tmpout_${ncfile} ${ncfile} done fi |
Wednesday, August 11, 2010
Arrays in Bash
| #!/bin/bash # # Usage: # FileNameMon.sh nmon # Copyright by http://scriptdemo.blogspot.com if [ $# -ne 1 ]; then sed -n '3,4p' FileNameMon.sh exit fi case $1 in 1) fnarray[0]="m01d05" fnarray[1]="m01d10" fnarray[2]="m01d15" fnarray[3]="m01d20" fnarray[4]="m01d25" fnarray[5]="m01d30" ;; 2) fnarray[0]="m02d04" fnarray[1]="m02d09" fnarray[2]="m02d14" fnarray[3]="m02d19" fnarray[4]="m02d24" fnarray[5]="m03d01" ;; 3) fnarray[0]="m03d06" fnarray[1]="m03d11" fnarray[2]="m03d16" fnarray[3]="m03d21" fnarray[4]="m03d26" fnarray[5]="m03d31" ;; 4) fnarray[0]="m04d05" fnarray[1]="m04d10" fnarray[2]="m04d15" fnarray[3]="m04d20" fnarray[4]="m04d25" fnarray[5]="m04d30" ;; 5) fnarray[0]="m05d05" fnarray[1]="m05d10" fnarray[2]="m05d15" fnarray[3]="m05d20" fnarray[4]="m05d25" fnarray[5]="m05d30" ;; 6) fnarray[0]="m06d04" fnarray[1]="m06d09" fnarray[2]="m06d14" fnarray[3]="m06d19" fnarray[4]="m06d24" fnarray[5]="m06d29" ;; 7) fnarray[0]="m07d04" fnarray[1]="m07d09" fnarray[2]="m07d14" fnarray[3]="m07d19" fnarray[4]="m07d24" fnarray[5]="m07d29" ;; 8) fnarray[0]="m08d03" fnarray[1]="m08d08" fnarray[2]="m08d13" fnarray[3]="m08d18" fnarray[4]="m08d23" fnarray[5]="m08d28" fnarray[6]="m09d02" ;; 9) fnarray[0]="m09d07" fnarray[1]="m09d12" fnarray[2]="m09d17" fnarray[3]="m09d22" fnarray[4]="m09d27" fnarray[5]="m10d02" ;; 10) fnarray[0]="m10d07" fnarray[1]="m10d12" fnarray[2]="m10d17" fnarray[3]="m10d22" fnarray[4]="m10d27" fnarray[5]="m11d01" ;; 11) fnarray[0]="m11d06" fnarray[1]="m11d11" fnarray[2]="m11d16" fnarray[3]="m11d21" fnarray[4]="m11d26" fnarray[5]="m12d01" ;; 12) fnarray[0]="m12d06" fnarray[1]="m12d11" fnarray[2]="m12d16" fnarray[3]="m12d21" fnarray[4]="m12d26" fnarray[5]="m12d31" ;; *) echo "Undefined month: $1" exit ;; esac echo ${fnarray[*]} |
mmddstr=(`FileNameMon.sh $month`) # will return an array with the month-day strings.
To get the array size, rather than the length of an array element, either ${#mmddstr[*]} or ${#mmddstr[@]} works well. One should notice that ${#mmddstr} shows something different but the max length of all the array elements.
Subscribe to:
Posts (Atom)