| 0,0 → 1,978 |
|---|
| PROGRAM z2s |
| |
| c Calculate secondary fields on z levels |
| c Michael Sprenger / Summer 2006 |
| |
| implicit none |
| |
| c --------------------------------------------------------------- |
| c Declaration of variables |
| c --------------------------------------------------------------- |
| |
| c Variables for output Z file : height level |
| character*80 cfn |
| real varmin(4),varmax(4),stag(4) |
| integer vardim(4) |
| real mdv |
| integer ndim |
| integer nx,ny,nz |
| real xmin,xmax,ymin,ymax,dx,dy |
| integer ntimes |
| real aklev(1000),bklev(1000) |
| real aklay(1000),bklay(1000) |
| real time |
| real pollon,pollat |
| integer idate(5) |
| integer nfields |
| character*80 field_nam(100) |
| real,allocatable, dimension (:,:,:,:) :: field_dat |
| real,allocatable, dimension (:,:,:) :: z3 |
| real,allocatable, dimension (:,:) :: x2,y2,f2 |
| real,allocatable, dimension (:,:,:) :: out |
| real,allocatable, dimension (:,:,:) :: inp |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| integer cdfid,cstid |
| |
| c Parameter file |
| character*80 fieldname |
| integer nfilt |
| character*80 ofn,gri |
| |
| c Auxiliary variables |
| integer ierr,stat |
| integer i,j,k,n |
| real,allocatable, dimension (:,:) :: out2 |
| character*80 vnam2(100) |
| integer nvars2 |
| |
| c --------------------------------------------------------------- |
| c Preparations |
| c --------------------------------------------------------------- |
| |
| print*,'*********************************************************' |
| print*,'* z2s *' |
| print*,'*********************************************************' |
| |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) fieldname |
| read(10,*) ofn |
| read(10,*) gri |
| read(10,*) nfilt |
| close(10) |
| |
| c Get grid description for Z file : height level |
| call cdfopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,'T',ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| nx =vardim(1) |
| ny =vardim(2) |
| nz =vardim(3) |
| xmin=varmin(1) |
| ymin=varmin(2) |
| call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr) |
| if (ierr.ne.0) goto 998 |
| call getgrid(cstid,dx,dy,ierr) |
| if (ierr.ne.0) goto 998 |
| xmax=xmin+real(nx-1)*dx |
| ymax=ymin+real(ny-1)*dy |
| call gettimes(cdfid,time,ntimes,ierr) |
| if (ierr.ne.0) goto 998 |
| call getstart(cstid,idate,ierr) |
| if (ierr.ne.0) goto 998 |
| call getpole(cstid,pollon,pollat,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Get a list of all variables on the GRID file |
| if (trim(gri).ne.trim(ofn)) then |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars2,vnam2,ierr) |
| if (ierr.ne.0) goto 998 |
| do i=1,nvars2 |
| vnam(nvars+i)=vnam2(i) |
| enddo |
| nvars=nvars+nvars2 |
| endif |
| |
| c Check whether calculation of <fieldname> is supported |
| if ( (fieldname.ne.'TH' ).and. |
| > (fieldname.ne.'NSQ') .and. |
| > (fieldname.ne.'PV' ) .and. |
| > (fieldname.ne.'RHO')) then |
| print*,'Variable not supported ',trim(fieldname) |
| stop |
| endif |
| |
| c Set dependencies |
| if (fieldname.eq.'TH') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| else if (fieldname.eq.'RHO') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| else if (fieldname.eq.'NSQ') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='Q' |
| else if (fieldname.eq.'PV') then |
| nfields=4 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| field_nam(3)='U' |
| field_nam(4)='V' |
| endif |
| |
| c Allocate memory |
| allocate(field_dat(nfields,nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating field_dat' |
| allocate(out(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating out' |
| allocate(inp(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating inp' |
| allocate(z3(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating z3' |
| allocate(x2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating x2' |
| allocate(y2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating y2' |
| allocate(f2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating f2' |
| |
| c Allocate auxiliary fields |
| allocate(out2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating out2' |
| |
| c Read X grid |
| varname='X' |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| print*,'Variable ',trim(varname),' missing... Stop' |
| stop |
| endif |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,x2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(gri) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Read Y grid |
| varname='Y' |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| print*,'Variable ',trim(varname),' missing... Stop' |
| stop |
| endif |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,y2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(gri) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Read Coriolis parametzer |
| varname='CORIOL' |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| print*,'Variable ',trim(varname),' missing... Stop' |
| stop |
| endif |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,f2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(gri) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Init height levels |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| z3(i,j,k)=aklay(k) |
| enddo |
| enddo |
| enddo |
| |
| c Load needed variables |
| do n=1,nfields |
| |
| c Check whether variable is available on file |
| varname=field_nam(n) |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| |
| c Variable is available on file |
| if (isok.eq.1) then |
| |
| call cdfopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,inp,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(ofn) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k)=inp(i,j,k) |
| enddo |
| enddo |
| enddo |
| |
| else |
| print*,'Variable missing : ',trim(varname) |
| stop |
| endif |
| |
| enddo |
| |
| c Change unit of pressure from hPa to Pa |
| n=0 |
| do i=1,nfields |
| if (field_nam(i).eq.'P') n=i |
| enddo |
| if (n.ne.0) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k)=100.*field_dat(n,i,j,k) |
| enddo |
| enddo |
| enddo |
| endif |
| |
| c ---------------------------------------------------------------- |
| c Calculations |
| c ---------------------------------------------------------------- |
| |
| c Call to the defining routines |
| if (fieldname.eq.'RHO') then |
| |
| print*,'C ',trim(fieldname) |
| call calc_rho (out, ! RHO |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! P |
| > nx,ny,nz,mdv) |
| |
| else if (fieldname.eq.'TH') then |
| |
| print*,'C ',trim(fieldname) |
| call calc_theta (out, ! TH |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! P |
| > nx,ny,nz,mdv) |
| |
| else if (fieldname.eq.'NSQ') then |
| |
| print*,'C ',trim(fieldname) |
| call calc_nsq (out, ! NSQ |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! Q |
| > z3, ! Z |
| > nx,ny,nz,mdv) |
| |
| else if (fieldname.eq.'PV') then |
| |
| print*,'C ',trim(fieldname) |
| call calc_pv (out, ! PV |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! P |
| > field_dat(3,:,:,:), ! U |
| > field_dat(4,:,:,:), ! V |
| > z3,x2,y2,f2, ! Z,X,Y,CORIOL |
| > nx,ny,nz,mdv) |
| |
| endif |
| |
| c Horizontal filtering of the resulting fields |
| print*,'F ',trim(fieldname) |
| do k=1,nz |
| |
| do i=1,nx |
| do j=1,ny |
| out2(i,j)=out(i,j,k) |
| enddo |
| enddo |
| do n=1,nfilt |
| call filt2d (out2,out2,nx,ny,1.,mdv,0,0,0,0) |
| enddo |
| |
| do i=1,nx |
| do j=1,ny |
| out(i,j,k)=out2(i,j) |
| enddo |
| enddo |
| |
| enddo |
| |
| c ---------------------------------------------------------------- |
| c Save result onto netcdf file |
| c ---------------------------------------------------------------- |
| |
| call cdfwopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname=fieldname |
| 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,out,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| |
| c ---------------------------------------------------------------- |
| c Exception handling |
| c ---------------------------------------------------------------- |
| |
| stop |
| |
| 998 print*,'Problem with netcdf file' |
| stop |
| |
| end |
| |
| |
| |
| c **************************************************************** |
| c * SUBROUTINE SECTION: AUXILIARY ROUTINES * |
| c **************************************************************** |
| |
| 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 |
| |
| |
| c **************************************************************** |
| c * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS * |
| c **************************************************************** |
| |
| c ----------------------------------------------------------------- |
| c Calculate density |
| c ----------------------------------------------------------------- |
| |
| subroutine calc_rho (rho3,t3,p3, |
| > nx,ny,nz,mdv) |
| |
| c Calculate the density <rho3> (in kg/m^3) from temperature <t3> |
| c (in deg C) and pressure <p3> (in hPa). |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real mdv |
| real rho3(nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real p3 (nx,ny,nz) |
| |
| c Declaration of physical constants |
| real eps |
| parameter (eps=0.01) |
| real rd |
| parameter (rd=287.) |
| real tzero |
| parameter (tzero=273.15) |
| |
| c Auxiliary variables |
| integer i,j,k |
| |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| |
| if ((abs(t3(i,j,k)-mdv).gt.eps).and. |
| > (abs(p3(i,j,k)-mdv).gt.eps)) then |
| |
| rho3(i,j,k)=1./rd*p3(i,j,k)/(t3(i,j,k)+tzero) |
| |
| else |
| rho3(i,j,k)=mdv |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| end |
| |
| |
| c ----------------------------------------------------------------- |
| c Calculate potential temperature |
| c ----------------------------------------------------------------- |
| |
| subroutine calc_theta (th3,t3,p3, |
| > nx,ny,nz,mdv) |
| |
| c Calculate potential temperature <th3> (in K) from temperature <t3> |
| c (in deg C) and pressure <p3> (in hPa). |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real mdv |
| real th3(nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real p3 (nx,ny,nz) |
| |
| c Declaration of physical constants |
| real rdcp,tzero,p0 |
| parameter (rdcp=0.286) |
| parameter (tzero=273.15) |
| parameter (p0=100000.) |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer i,j,k |
| |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| |
| if ((abs(t3(i,j,k)-mdv).gt.eps).and. |
| > (abs(p3(i,j,k)-mdv).gt.eps)) then |
| |
| th3(i,j,k)=(t3(i,j,k)+tzero)*( (p0/p3(i,j,k))**rdcp ) |
| |
| else |
| th3(i,j,k)=mdv |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| end |
| |
| c ----------------------------------------------------------------- |
| c Calculate stratification |
| c ----------------------------------------------------------------- |
| |
| subroutine calc_nsq (nsq3,t3,q3, |
| > z3,nx,ny,nz,mdv) |
| |
| c Calculate stratification <nsq3> on the model grid. The grid is |
| c specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number |
| c of vertical levels is <nz>. The input field are: temperature <t3>, |
| c specific humidity <q3>, height <z3> and horizontal grid <x2>, <y2>. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real nsq3(nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real q3 (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real mdv |
| |
| c Physical and numerical parameters |
| real scale |
| parameter (scale=1.) |
| real g |
| parameter (g=9.80665) |
| real eps |
| parameter (eps=0.01) |
| real tzero |
| parameter (tzero=273.15) |
| real kappa |
| parameter (kappa=0.6078) |
| real cp |
| parameter (cp=1005.) |
| real zerodiv |
| parameter (zerodiv=0.0000000001) |
| real R |
| parameter (R=287.) |
| |
| c Auxiliary variables |
| real tv (nx,ny,nz) |
| real dtvdz(nx,ny,nz) |
| integer i,j,k |
| real scaledz |
| |
| c Calculate 3d virtual temperature |
| do k=1,nz |
| do i=1,nx |
| do j=1,ny |
| if ((abs(t3(i,j,k)-mdv).gt.eps).and. |
| > (abs(q3(i,j,k)-mdv).gt.eps)) |
| > then |
| tv(i,j,k) = (t3(i,j,k)+tzero)*(1.+kappa*q3(i,j,k)) |
| else |
| tv(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c Vertical derivative of virtual temperature |
| call deriv (dtvdz,tv,'z',z3,x2,y2,nx,ny,nz,mdv) |
| |
| c Stratification |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(dtvdz(i,j,k)-mdv).gt.eps).and. |
| > (abs(tv (i,j,k)-mdv).gt.eps)) |
| > then |
| nsq3(i,j,k) = g/tv(i,j,k) * (dtvdz(i,j,k) + g/cp) |
| else |
| nsq3(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| |
| end |
| |
| c ----------------------------------------------------------------- |
| c Calculate potential vorticity |
| c ----------------------------------------------------------------- |
| |
| subroutine calc_pv (pv3,t3,p3,u3,v3, |
| > z3,x2,y2,f2,nx,ny,nz,mdv) |
| |
| c Calculate potential vorticity <pv3> on the model grid. The grid is |
| c specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number |
| c of vertical levels is <nz>. The input field are: potential temperature |
| c <th3>, horizontal wind <u3> and <v3>, density <rho3>, height <z3>. |
| c The terms including the vertical velocity are neglected in the calculation |
| c of the PV. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real pv3 (nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real p3 (nx,ny,nz) |
| real u3 (nx,ny,nz) |
| real v3 (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real f2 (nx,ny) |
| |
| real mdv |
| |
| c Physical and numerical parameters |
| real scale |
| parameter (scale=1.E6) |
| real omega |
| parameter (omega=7.292E-5) |
| real pi180 |
| parameter (pi180=3.141592654/180.) |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| real dtdz(nx,ny,nz) |
| real dtdx(nx,ny,nz) |
| real dtdy(nx,ny,nz) |
| real dudz(nx,ny,nz) |
| real dvdz(nx,ny,nz) |
| real dvdx(nx,ny,nz) |
| real dudy(nx,ny,nz) |
| real rho3(nx,ny,nz) |
| real th3 (nx,ny,nz) |
| |
| integer i,j,k |
| |
| c Calculate density and potential temperature |
| call calc_rho (rho3,t3,p3,nx,ny,nz,mdv) |
| call calc_theta (th3 ,t3,p3,nx,ny,nz,mdv) |
| |
| c Calculate needed derivatives |
| call deriv (dudz, u3,'z',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dvdz, v3,'z',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dtdz,th3,'z',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv) |
| |
| c Calculate potential vorticity |
| do j=1,ny |
| do i=1,nx |
| do k=1,nz |
| |
| c Evaluate PV formula with missing data check |
| if ((abs(dtdz(i,j,k)-mdv).gt.eps).and. |
| > (abs(dudz(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdy(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdz(i,j,k)-mdv).gt.eps).and. |
| > (abs(rho3(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dudy(i,j,k)-mdv).gt.eps)) then |
| |
| pv3(i,j,k)=scale*1./rho3(i,j,k)* |
| > ((dvdx(i,j,k)-dudy(i,j,k)+f2(i,j))*dtdz(i,j,k) |
| > +dudz(i,j,k)*dtdy(i,j,k) |
| > -dvdz(i,j,k)*dtdx(i,j,k)) |
| else |
| pv3(i,j,k)=mdv |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| end |
| |
| |
| c **************************************************************** |
| c * SUBROUTINE SECTION: GRID HANDLING * |
| c **************************************************************** |
| |
| c ----------------------------------------------------------------- |
| c Horizontal and vertical derivatives for 3d fields |
| c ----------------------------------------------------------------- |
| |
| subroutine deriv (df,f,direction, |
| > z3,x2,y2,nx,ny,nz,mdv) |
| |
| c Calculate horizontal and vertical derivatives of the 3d field <f>. |
| c The direction of the derivative is specified in <direction> |
| c 'x','y' : Horizontal derivative in x and y direction |
| c 'p','z','t','m' : Vertical derivative (pressure, height, theta, model) |
| c The 3d field <z3> specifies the isosurfaces along which the horizontal |
| c derivatives are calculated or the levels for the vertical derivatives. |
| |
| implicit none |
| |
| c Input and output parameters |
| integer nx,ny,nz |
| real df (nx,ny,nz) |
| real f (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| character direction |
| real mdv |
| |
| c Numerical and physical parameters |
| real pi180 |
| parameter (pi180=3.141592654/180.) |
| real zerodiv |
| parameter (zerodiv=0.00000001) |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer i,j,k |
| real vmin,vmax |
| real scale,lat |
| real vu,vl,vuvl,vlvu |
| integer o,u,w,e,n,s |
| |
| c Vertical derivative |
| if ((direction.eq.'z').or. |
| > (direction.eq.'th').or. |
| > (direction.eq.'p').or. |
| > (direction.eq.'m').and. |
| > (nz.gt.1)) then |
| |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| |
| o=k+1 |
| if (o.gt.nz) o=nz |
| u=k-1 |
| if (u.lt.1) u=1 |
| |
| if ((abs(f(i,j,o)-mdv).gt.eps).and. |
| > (abs(f(i,j,u)-mdv).gt.eps).and. |
| > (abs(f(i,j,k)-mdv).gt.eps)) then |
| |
| vu = z3(i,j,k)-z3(i,j,o) |
| vl = z3(i,j,u)-z3(i,j,k) |
| vuvl = vu/(vl+zerodiv) |
| vlvu = 1./(vuvl+zerodiv) |
| |
| df(i,j,k) = 1./(vu+vl) |
| > * (vuvl*(f(i,j,u)-f(i,j,k)) |
| > + vlvu*(f(i,j,k)-f(i,j,o))) |
| |
| else |
| df(i,j,k) = mdv |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| c Horizontal derivative in the y direction: 3d |
| elseif (direction.eq.'y') then |
| |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| |
| s=j-1 |
| if (s.lt.1) s=1 |
| n=j+1 |
| if (n.gt.ny) n=ny |
| |
| if ((abs(f(i,n,k)-mdv).gt.eps).and. |
| > (abs(f(i,j,k)-mdv).gt.eps).and. |
| > (abs(f(i,s,k)-mdv).gt.eps)) then |
| |
| vu = 1000.*(y2(i,j)-y2(i,n)) |
| vl = 1000.*(y2(i,s)-y2(i,j)) |
| vuvl = vu/(vl+zerodiv) |
| vlvu = 1./(vuvl+zerodiv) |
| |
| df(i,j,k) = 1./(vu+vl) |
| > * (vuvl*(f(i,s,k)-f(i,j,k)) |
| > + vlvu*(f(i,j,k)-f(i,n,k))) |
| |
| else |
| df(i,j,k) = mdv |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| c Horizontal derivative in the x direction: 3d |
| elseif (direction.eq.'x') then |
| |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| |
| w=i-1 |
| if (w.lt.1) w=1 |
| e=i+1 |
| if (e.gt.nx) e=nx |
| |
| if ((abs(f(w,j,k)-mdv).gt.eps).and. |
| > (abs(f(i,j,k)-mdv).gt.eps).and. |
| > (abs(f(e,j,k)-mdv).gt.eps)) then |
| |
| vu = 1000.*(x2(i,j)-x2(e,j)) |
| vl = 1000.*(x2(w,j)-x2(i,j)) |
| vuvl = vu/(vl+zerodiv) |
| vlvu = 1./(vuvl+zerodiv) |
| |
| df(i,j,k) = 1./(vu+vl) |
| > * (vuvl*(f(w,j,k)-f(i,j,k)) |
| > + vlvu*(f(i,j,k)-f(e,j,k))) |
| |
| else |
| df(i,j,k) = mdv |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| c Undefined direction for derivative |
| else |
| |
| print*,'Invalid direction of derivative... Stop' |
| stop |
| |
| endif |
| |
| end |
| |
| c ----------------------------------------------------------------- |
| c Horizontal filter |
| c ----------------------------------------------------------------- |
| |
| subroutine filt2d (a,af,nx,ny,fil,misdat, |
| & iperx,ipery,ispol,inpol) |
| |
| c Apply a conservative diffusion operator onto the 2d field a, |
| c with full missing data checking. |
| c |
| c a real inp array to be filtered, dimensioned (nx,ny) |
| c af real out filtered array, dimensioned (nx,ny), can be |
| c equivalenced with array a in the calling routine |
| c f1 real workarray, dimensioned (nx+1,ny) |
| c f2 real workarray, dimensioned (nx,ny+1) |
| c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1 |
| c corresponds to one application of the linear filter. |
| c misdat real inp missing-data value, a(i,j)=misdat indicates that |
| c the corresponding value is not available. The |
| c misdat-checking can be switched off with with misdat=0. |
| c iperx int inp periodic boundaries in the x-direction (1=yes,0=no) |
| c ipery int inp periodic boundaries in the y-direction (1=yes,0=no) |
| c inpol int inp northpole at j=ny (1=yes,0=no) |
| c ispol int inp southpole at j=1 (1=yes,0=no) |
| c |
| c Christoph Schaer, 1993 |
| |
| c argument declaration |
| integer nx,ny |
| real a(nx,ny),af(nx,ny),fil,misdat |
| integer iperx,ipery,inpol,ispol |
| |
| c local variable declaration |
| integer i,j,is |
| real fh |
| real f1(nx+1,ny),f2(nx,ny+1) |
| |
| c compute constant fh |
| fh=0.125*fil |
| |
| c compute fluxes in x-direction |
| if (misdat.eq.0.) then |
| do j=1,ny |
| do i=2,nx |
| f1(i,j)=a(i-1,j)-a(i,j) |
| enddo |
| enddo |
| else |
| do j=1,ny |
| do i=2,nx |
| if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then |
| f1(i,j)=0. |
| else |
| f1(i,j)=a(i-1,j)-a(i,j) |
| endif |
| enddo |
| enddo |
| endif |
| if (iperx.eq.1) then |
| c do periodic boundaries in the x-direction |
| do j=1,ny |
| f1(1,j)=f1(nx,j) |
| f1(nx+1,j)=f1(2,j) |
| enddo |
| else |
| c set boundary-fluxes to zero |
| do j=1,ny |
| f1(1,j)=0. |
| f1(nx+1,j)=0. |
| enddo |
| endif |
| |
| c compute fluxes in y-direction |
| if (misdat.eq.0.) then |
| do j=2,ny |
| do i=1,nx |
| f2(i,j)=a(i,j-1)-a(i,j) |
| enddo |
| enddo |
| else |
| do j=2,ny |
| do i=1,nx |
| if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then |
| f2(i,j)=0. |
| else |
| f2(i,j)=a(i,j-1)-a(i,j) |
| endif |
| enddo |
| enddo |
| endif |
| c set boundary-fluxes to zero |
| do i=1,nx |
| f2(i,1)=0. |
| f2(i,ny+1)=0. |
| enddo |
| if (ipery.eq.1) then |
| c do periodic boundaries in the x-direction |
| do i=1,nx |
| f2(i,1)=f2(i,ny) |
| f2(i,ny+1)=f2(i,2) |
| enddo |
| endif |
| if (iperx.eq.1) then |
| if (ispol.eq.1) then |
| c do south-pole |
| is=(nx-1)/2 |
| do i=1,nx |
| f2(i,1)=-f2(mod(i-1+is,nx)+1,2) |
| enddo |
| endif |
| if (inpol.eq.1) then |
| c do north-pole |
| is=(nx-1)/2 |
| do i=1,nx |
| f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny) |
| enddo |
| endif |
| endif |
| |
| c compute flux-convergence -> filter |
| if (misdat.eq.0.) then |
| do j=1,ny |
| do i=1,nx |
| af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1)) |
| enddo |
| enddo |
| else |
| do j=1,ny |
| do i=1,nx |
| if (a(i,j).eq.misdat) then |
| af(i,j)=misdat |
| else |
| af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1)) |
| endif |
| enddo |
| enddo |
| endif |
| end |
| Property changes: |
|---|
| Added: svn:executable |
|---|