| 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 |
|---|