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 |