Subversion Repositories pvinversion.ecmwf

Compare Revisions

No changes between revisions

Ignore whitespace Rev 2 → Rev 3

/trunk/prep/coastline.f
0,0 → 1,528
PROGRAM coastline
c *********************************************************************
c * Add the coastline and land sea mask to the REF file *
c * Michael Sprenger / Winter 2006,07 *
c *********************************************************************
 
implicit none
 
c --------------------------------------------------------------------
c Declaration of parameters and variables
c --------------------------------------------------------------------
 
c Input parameters
character*80 coastfile
character*80 reffile
real clon,clat,crot
c Output arrays
integer nmax
parameter (nmax=100000)
real,allocatable, dimension (:,:) :: landsea
real,allocatable, dimension (:,:) :: xcor,ycor
real xpos(nmax),ypos(nmax)
real lon(nmax), lat(nmax)
real rlon(nmax),rlat(nmax)
integer ncoast
integer rnx,rny
real rxmin,rymin,rxmax,rymax
real rdx,rdy
 
c Auxiliary variables
integer i,j
character*80 name
integer cdfid,cstid
character*80 cfn
integer ierr,stat
character*80 varname
real varmin(4),varmax(4),stag(4)
integer vardim(4),ndim
real mdv
integer isok
integer nvars
integer ntimes
real time
character*80 vnam(100)
 
c --------------------------------------------------------------------
c Preparations
c --------------------------------------------------------------------
 
print*,'*********************************************************'
print*,'* coastline *'
print*,'*********************************************************'
 
 
c Read parameter file
open(10,file='fort.10')
read(10,*) reffile
read(10,*) coastfile
read(10,*) name,clon
if ( trim(name).ne.'CLON' ) stop
read(10,*) name,clat
if ( trim(name).ne.'CLAT' ) stop
read(10,*) name,crot
if ( trim(name).ne.'CROT' ) stop
close(10)
 
print*
print*,trim(reffile)
print*,trim(coastfile)
print*,clon,clat,crot
 
c Read grid parameters
call cdfopn(reffile,cdfid,ierr)
if (ierr.ne.0) goto 998
call getvars(cdfid,nvars,vnam,ierr)
if (ierr.ne.0) goto 998
varname='X'
call getdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
call gettimes(cdfid,time,ntimes,ierr)
if (ierr.ne.0) goto 998
call clscdf(cdfid,ierr)
rnx=vardim(1)
rny=vardim(2)
rxmin=varmin(1)
rymin=varmin(2)
rxmax=varmax(1)
rymax=varmax(2)
rdx=(rxmax-rxmin)/real(rnx-1)
rdy=(rymax-rymin)/real(rny-1)
 
c Allocate memory
allocate(landsea(rnx,rny),stat=stat)
if (stat.ne.0) print*,'*** error allocating array landsea ***'
allocate(xcor(rnx,rny),stat=stat)
if (stat.ne.0) print*,'*** error allocating array xcor ***'
allocate(ycor(rnx,rny),stat=stat)
if (stat.ne.0) print*,'*** error allocating array ycor ***'
 
c Read data
call cdfopn(reffile,cdfid,ierr)
if (ierr.ne.0) goto 998
varname='X'
call getdat(cdfid,varname,time,1,xcor,ierr)
if (ierr.ne.0) goto 998
varname='Y'
call getdat(cdfid,varname,time,1,ycor,ierr)
if (ierr.ne.0) goto 998
call clscdf(cdfid,ierr)
 
c Read the coastfile
open(10,file=coastfile)
ncoast=1
100 read(10,*,end=200) lon(ncoast),lat(ncoast)
ncoast=ncoast+1
goto 100
 
200 close(10)
ncoast=ncoast-1
 
c Handle missing data values
mdv=-100000.
do i=ncoast,nmax
lon(i)=mdv
lat(i)=mdv
enddo
 
c --------------------------------------------------------------------
c Transform coast(lat/lon) to coast(x/y)
c --------------------------------------------------------------------
call getenvir (clon,clat,crot,
> xcor,ycor,rnx,rny,rxmin,rymin,rdx,rdy,mdv,
> lon,lat,rlon,rlat,xpos,ypos,nmax)
 
c --------------------------------------------------------------------
c Write output
c --------------------------------------------------------------------
 
c Open output file and set some general parameters
call cdfwopn(reffile,cdfid,ierr)
if (ierr.ne.0) goto 998
vardim(1)=1
vardim(2)=1
vardim(3)=nmax
ndim=3
 
c Write LAT
isok=0
varname='COAST_LAT'
call check_varok(isok,varname,vnam,nvars)
if (isok.eq.0) then
call putdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
endif
call putdat(cdfid,varname,time,0,lat,ierr)
if (ierr.ne.0) goto 998
print*,'W ',trim(varname),' ',trim(reffile)
 
c Write LON
isok=0
varname='COAST_LON'
call check_varok(isok,varname,vnam,nvars)
if (isok.eq.0) then
call putdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
endif
call putdat(cdfid,varname,time,0,lon,ierr)
if (ierr.ne.0) goto 998
print*,'W ',trim(varname),' ',trim(reffile)
 
c Write RLAT
isok=0
varname='COAST_RLAT'
call check_varok(isok,varname,vnam,nvars)
if (isok.eq.0) then
call putdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
endif
call putdat(cdfid,varname,time,0,rlat,ierr)
if (ierr.ne.0) goto 998
print*,'W ',trim(varname),' ',trim(reffile)
 
c Write RLON
isok=0
varname='COAST_RLON'
call check_varok(isok,varname,vnam,nvars)
if (isok.eq.0) then
call putdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
endif
call putdat(cdfid,varname,time,0,rlon,ierr)
if (ierr.ne.0) goto 998
print*,'W ',trim(varname),' ',trim(reffile)
 
c Write Y
isok=0
varname='COAST_Y'
call check_varok(isok,varname,vnam,nvars)
if (isok.eq.0) then
call putdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
endif
call putdat(cdfid,varname,time,0,ypos,ierr)
if (ierr.ne.0) goto 998
print*,'W ',trim(varname),' ',trim(reffile)
 
c Write X
isok=0
varname='COAST_X'
call check_varok(isok,varname,vnam,nvars)
if (isok.eq.0) then
call putdef(cdfid,varname,ndim,mdv,vardim,
> varmin,varmax,stag,ierr)
if (ierr.ne.0) goto 998
endif
call putdat(cdfid,varname,time,0,xpos,ierr)
if (ierr.ne.0) goto 998
print*,'W ',trim(varname),' ',trim(reffile)
 
c Close output file
call clscdf(cdfid,ierr)
if (ierr.ne.0) goto 998
 
 
c --------------------------------------------------------------------
c Exception handling
c --------------------------------------------------------------------
stop
 
998 print*,'Problems with input netcdf file'
stop
 
END
 
c --------------------------------------------------------------------------------
c Subroutine to get environment of strcof
c --------------------------------------------------------------------------------
 
SUBROUTINE getenvir (clon,clat,rotation,
> xcor,ycor,rnx,rny,rxmin,rymin,rdx,rdy,mdv,
> lon,lat,rlon,rlat,xpos,ypos,n)
 
c Rotate from a local quasi-cartesian coordiante system into lat/lon coordinate
c system.
 
implicit none
 
c Declaration of input and output parameters
integer n
real clon,clat,rotation
real lon(n), lat(n)
real rlon(n),rlat(n)
real xpos(n),ypos(n)
integer rnx,rny
real xcor(rnx,rny),ycor(rnx,rny)
real rxmin,rymin,rdx,rdy
real mdv
 
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 Auxiliary variables
real pollon,pollat
integer i,j
real xind,yind
real rindx,rindy
integer indx,indy,indr,indu
real frac0i,frac0j,frac1i,frac1j
real rlon1(n),rlat1(n)
 
c Externals
real lmtolms,phtophs
external lmtolms,phtophs
 
c First rotation
pollon=clon-180.
if (pollon.lt.-180.) pollon=pollon+360.
pollat=90.-clat
do i=1,n
 
if ( (abs(lon(i)-mdv).gt.eps).and.
> (abs(lat(i)-mdv).gt.eps) ) then
 
c First rotation
pollon=clon-180.
if (pollon.lt.-180.) pollon=pollon+360.
pollat=90.-clat
rlon1(i)=lmtolms(lat(i),lon(i),pollat,pollon)
rlat1(i)=phtophs(lat(i),lon(i),pollat,pollon)
 
c Second coordinate transformation
pollon=-180.
pollat=90.+rotation
rlon(i)=90.+lmtolms(rlat1(i),rlon1(i)-90.,pollat,pollon)
rlat(i)=phtophs(rlat1(i),rlon1(i)-90.,pollat,pollon)
 
c Get rotated latitude and longitude
100 if (rlon(i).lt.rxmin) then
rlon(i)=rlon(i)+360.
goto 100
endif
102 if (rlon(i).gt.(rxmin+real(rnx-1)*rdx)) then
rlon(i)=rlon(i)-360.
goto 102
endif
if ( (rlon(i).lt.rxmin).or.
> (rlon(i).gt.(rxmin+real(rnx-1)*rdx)).or.
> (rlat(i).lt.rymin).or.
> (rlat(i).gt.(rymin+real(rny-1)*rdy)) ) then
rlon(i)=mdv
rlat(i)=mdv
endif
 
c Get x and y coordinates
rindx=(rlon(i)-rxmin)/rdx+1.
rindy=(rlat(i)-rymin)/rdy+1.
indx=int(rindx)
indy=int(rindy)
if ((indx.gt.1).and.(indx.lt.rnx).and.
> (indy.gt.1).and.(indy.lt.rny)) then
 
indr=indx+1
if (indr.gt.rnx) indr=1
indu=indy+1
if (indu.gt.rny) indu=rny
 
frac0i=rindx-float(indx)
frac0j=rindy-float(indy)
frac1i=1.-frac0i
frac1j=1.-frac0j
xpos(i) = xcor(indx ,indy) * frac1i * frac1j
& + xcor(indx ,indu) * frac1i * frac0j
& + xcor(indr ,indy) * frac0i * frac1j
& + xcor(indr ,indu) * frac0i * frac0j
ypos(i) = ycor(indx ,indy) * frac1i * frac1j
& + ycor(indx ,indu) * frac1i * frac0j
& + ycor(indr ,indy) * frac0i * frac1j
& + ycor(indr ,indu) * frac0i * frac0j
 
else
 
xpos(i)=mdv
ypos(i)=mdv
 
endif
 
else
rlon(i)=mdv
rlat(i)=mdv
xpos(i)=mdv
ypos(i)=mdv
endif
 
enddo
 
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 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