| 0,0 → 1,1567 |
|---|
| PROGRAM add2p |
| |
| c ********************************************************************* |
| c * Add output from the PV inversion to the original P file * |
| c * Zonal wind (U), meridional wind (V) and temperature (T) are * |
| c * modified. All other fields are not affected by the programme * |
| c * * |
| c * Michael Sprenger / Summer, Autumn 2006 |
| c ********************************************************************* |
| |
| implicit none |
| |
| c --------------------------------------------------------------------- |
| c Declaration of parameters and variables |
| c --------------------------------------------------------------------- |
| |
| c Input parameters |
| character*80 ml_fn |
| character*80 zm_fn |
| character*80 or_fn |
| real transxy |
| |
| c Variables for input P file : model level |
| real ml_varmin(4),ml_varmax(4),ml_stag(4) |
| integer ml_vardim(4) |
| real ml_mdv |
| integer ml_ndim |
| integer ml_nx,ml_ny,ml_nz |
| real ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy |
| integer ml_ntimes |
| real ml_aklev(500),ml_bklev(500) |
| real ml_aklay(500),ml_bklay(500) |
| real ml_time |
| real ml_pollon,ml_pollat |
| integer ml_idate(5) |
| integer ml_nvars |
| character*80 ml_vnam(80) |
| real,allocatable, dimension (:,:) :: ml_ps,ml_oro |
| real,allocatable, dimension (:,:,:) :: ml_t3,ml_u3,ml_v3,ml_z3 |
| real,allocatable, dimension (:,:,:) :: ml_p3,ml_tv3,ml_q3,ml_p0 |
| integer,allocatable, dimension (:,:,:) :: flag_ml |
| |
| c Variables for GEO.MOD |
| real zm_varmin(4),zm_varmax(4),zm_stag(4) |
| integer zm_vardim(4) |
| real zm_mdv |
| integer zm_ndim |
| integer zm_nx,zm_ny,zm_nz |
| real zm_xmin,zm_xmax,zm_ymin,zm_ymax,zm_dx,zm_dy |
| integer zm_ntimes |
| real zm_aklev(500),zm_bklev(500) |
| real zm_aklay(500),zm_bklay(500) |
| real zm_time |
| real zm_pollon,zm_pollat |
| integer zm_idate(5) |
| integer zm_nvars |
| character*80 zm_vnam(80) |
| real,allocatable, dimension (:,:,:) :: zm_u3,zm_v3,zm_t3,zm_z3 |
| real,allocatable, dimension (:,:,:) :: zm_p3 |
| real,allocatable, dimension (:,:) :: zm_rlat,zm_rlon |
| integer,allocatable, dimension (:,:,:) :: flag_zm |
| |
| c Variables for GEO.ORG |
| real or_varmin(4),or_varmax(4),or_stag(4) |
| integer or_vardim(4) |
| real or_mdv |
| integer or_ndim |
| integer or_nx,or_ny,or_nz |
| real or_xmin,or_xmax,or_ymin,or_ymax,or_dx,or_dy |
| integer or_ntimes |
| real or_time |
| real or_pollon,or_pollat |
| integer or_idate(5) |
| integer or_nvars |
| character*80 or_vnam(80) |
| real,allocatable, dimension (:,:,:) :: or_p3 |
| real,allocatable, dimension (:,:,:) :: or_u3,or_v3,or_t3 |
| |
| c Anomalies |
| real,allocatable, dimension (:,:,:) :: an_p3,an_u3,an_v3,an_t3 |
| |
| c Array with the weighting function |
| real,allocatable, dimension (:,:) :: weight,dist |
| |
| c Flag for test mode |
| integer test |
| parameter (test=1) |
| |
| c Flag for Poisson filling |
| integer poisson |
| parameter (poisson=1) |
| |
| c Physical and numerical parameters |
| real eps ! Numerical epsilon |
| parameter (eps=0.001) |
| real dpextra ! Allowed range for extrapolation (% of pressure) |
| parameter (dpextra=0.1) |
| real g ! Earth's gravity |
| parameter (g=9.80665) |
| real tzero ! Conversion Celius/Kelvin |
| parameter (tzero=273.15) |
| real kappa ! Kappa (Definition of potential temperature) |
| parameter (kappa=0.6078) |
| real zerodiv ! Zero division handling |
| parameter (zerodiv=0.0000000001) |
| real dpmin ! Pressure step for hydrostatic equation |
| parameter (dpmin=10.) |
| real rdg ! Ratio gas constant / Earth's gravity |
| parameter (rdg=29.271) |
| |
| c Auxiliray variables |
| integer i,j,k,l |
| character*80 varname |
| integer isok |
| integer stat |
| integer cdfid,ierr,cstid |
| character*80 cfn |
| real p1(1000),z1(1000),f1(1000),tv1(1000) |
| real spline_f1(1000),spline_tv1(1000) |
| integer n1 |
| real hh |
| real pmin,pmax |
| real boundlat(100000),boundlon(1000000) |
| integer nbound |
| integer il,ir,ju,jd,im,jm |
| real lat,lon |
| real distmin,distpos |
| character*80 name |
| real pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ff |
| integer lmin,n |
| real mean,nmean |
| real psfc |
| integer count0,count1 |
| |
| c Externals |
| real sdis,kink |
| external sdis,kink |
| |
| |
| c ----------------------------------------------------------------- |
| c Read input fields |
| c ----------------------------------------------------------------- |
| |
| print*,'*********************************************************' |
| print*,'* add2p *' |
| print*,'*********************************************************' |
| print* |
| |
| c Read in the parameter file |
| open(10,file='fort.10') |
| read(10,*) ml_fn |
| read(10,*) zm_fn |
| read(10,*) or_fn |
| read(10,*) name,transxy |
| if ( trim(name).ne.'TRANS_XY' ) stop |
| read(10,*) name,psfc |
| if ( trim(name).ne.'PS_CHANGE' ) stop |
| close(10) |
| print*,trim(ml_fn) |
| print*,trim(zm_fn) |
| print*,trim(or_fn) |
| print* |
| |
| c Get grid description for P file : model level |
| call cdfopn(ml_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,ml_nvars,ml_vnam,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 |
| isok=0 |
| varname='T' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.ne.1) goto 998 |
| call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| ml_nx =ml_vardim(1) |
| ml_ny =ml_vardim(2) |
| ml_nz =ml_vardim(3) |
| ml_xmin=ml_varmin(1) |
| ml_ymin=ml_varmin(2) |
| call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr) |
| if (ierr.ne.0) goto 998 |
| call getgrid(cstid,ml_dx,ml_dy,ierr) |
| ml_xmax=ml_xmin+real(ml_nx-1)*ml_dx |
| ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy |
| if (ierr.ne.0) goto 998 |
| call gettimes(cdfid,ml_time,ml_ntimes,ierr) |
| if (ierr.ne.0) goto 998 |
| call getstart(cstid,ml_idate,ierr) |
| if (ierr.ne.0) goto 998 |
| call getpole(cstid,ml_pollon,ml_pollat,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 grid description for MOD file |
| call cdfopn(zm_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| call getvars(cdfid,zm_nvars,zm_vnam,ierr) |
| if (ierr.ne.0) goto 997 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 997 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 997 |
| isok=0 |
| varname='T' |
| call check_varok(isok,varname,zm_vnam,zm_nvars) |
| if (isok.ne.1) goto 997 |
| call getdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim, |
| > zm_varmin,zm_varmax,zm_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| zm_nx =zm_vardim(1) |
| zm_ny =zm_vardim(2) |
| zm_nz =zm_vardim(3) |
| zm_xmin=zm_varmin(1) |
| zm_ymin=zm_varmin(2) |
| call getlevs(cstid,zm_nz,zm_aklev,zm_bklev,zm_aklay,zm_bklay,ierr) |
| if (ierr.ne.0) goto 997 |
| call getgrid(cstid,zm_dx,zm_dy,ierr) |
| if (ierr.ne.0) goto 997 |
| zm_xmax=zm_xmin+real(zm_nx-1)*zm_dx |
| zm_ymax=zm_ymin+real(zm_ny-1)*zm_dy |
| call gettimes(cdfid,zm_time,zm_ntimes,ierr) |
| if (ierr.ne.0) goto 997 |
| call getstart(cstid,zm_idate,ierr) |
| if (ierr.ne.0) goto 997 |
| call getpole(cstid,zm_pollon,zm_pollat,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Get grid description for ORG file |
| call cdfopn(or_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| call getvars(cdfid,or_nvars,or_vnam,ierr) |
| if (ierr.ne.0) goto 997 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 997 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 997 |
| isok=0 |
| varname='P' |
| call check_varok(isok,varname,or_vnam,or_nvars) |
| if (isok.ne.1) goto 997 |
| call getdef(cdfid,varname,or_ndim,or_mdv,or_vardim, |
| > or_varmin,or_varmax,or_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| or_nx =or_vardim(1) |
| or_ny =or_vardim(2) |
| or_nz =or_vardim(3) |
| or_xmin=or_varmin(1) |
| or_ymin=or_varmin(2) |
| call getgrid(cstid,or_dx,or_dy,ierr) |
| if (ierr.ne.0) goto 997 |
| or_xmax=or_xmin+real(or_nx-1)*or_dx |
| or_ymax=or_ymin+real(or_ny-1)*or_dy |
| call gettimes(cdfid,or_time,or_ntimes,ierr) |
| if (ierr.ne.0) goto 997 |
| call getstart(cstid,or_idate,ierr) |
| if (ierr.ne.0) goto 997 |
| call getpole(cstid,or_pollon,or_pollat,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Consistency check for the grids |
| if ( (abs(ml_xmin-zm_xmin).gt.eps).or. |
| > (abs(ml_xmax-zm_xmax).gt.eps).or. |
| > (abs(ml_ymin-zm_ymin).gt.eps).or. |
| > (abs(ml_ymax-zm_ymax).gt.eps).or. |
| > (abs(ml_dx -zm_dx ).gt.eps).or. |
| > (abs(ml_dy -zm_dy ).gt.eps).or. |
| > (abs(ml_time-zm_time).gt.eps).or. |
| > (abs(ml_xmin-or_xmin).gt.eps).or. |
| > (abs(ml_xmax-or_xmax).gt.eps).or. |
| > (abs(ml_ymin-or_ymin).gt.eps).or. |
| > (abs(ml_ymax-or_ymax).gt.eps).or. |
| > (abs(ml_dx -or_dx ).gt.eps).or. |
| > (abs(ml_dy -or_dy ).gt.eps) ) |
| >then |
| print*,'P, GEO.MOD, GEO.ORG grids not consistent... Stop' |
| print*,ml_xmin,zm_xmin,or_xmin |
| print*,ml_ymin,zm_ymin,or_ymin |
| print*,ml_dx, zm_dx ,or_dx |
| print*,ml_dy, zm_dy ,or_dy |
| print*,ml_time,zm_time,or_time |
| stop |
| endif |
| |
| c Allocate memory for input and output P file |
| allocate(ml_p3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_p3' |
| allocate(ml_u3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_u3' |
| allocate(ml_v3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_v3' |
| allocate(ml_t3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_t3' |
| allocate(ml_tv3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_tv3' |
| allocate(ml_z3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_z3' |
| allocate(ml_q3 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_q3' |
| allocate(ml_ps (ml_nx,ml_ny ),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_ps' |
| allocate(ml_oro (ml_nx,ml_ny ),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_oro' |
| allocate(ml_p0 (ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating ml_p0' |
| |
| c Allocate memory for input GEO.MOD file |
| allocate(zm_p3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_p3' |
| allocate(zm_u3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_u3' |
| allocate(zm_v3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_v3' |
| allocate(zm_t3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_t3' |
| allocate(zm_z3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_z3' |
| allocate(zm_rlat(zm_nx,zm_ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_rlat' |
| allocate(zm_rlon(zm_nx,zm_ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_rlon' |
| |
| c Allocate memory for input GEO.ORG file |
| allocate(or_p3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_p3' |
| allocate(or_t3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_t3' |
| allocate(or_u3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_u3' |
| allocate(or_v3(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating zm_v3' |
| |
| c Allocate memory for weighting function, flag and anomalies |
| allocate(weight(zm_nx,zm_ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating weight' |
| allocate(dist (zm_nx,zm_ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating dist' |
| allocate(flag_zm(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating flag' |
| allocate(flag_ml(zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating flag' |
| allocate(an_p3 (zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating an_p3' |
| allocate(an_t3 (zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating an_t3' |
| allocate(an_u3 (zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating an_u3' |
| allocate(an_v3 (zm_nx,zm_ny,zm_nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating an_v3' |
| |
| c Read variables from P file : model levels |
| call cdfopn(ml_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| varname='T' |
| call getdat(cdfid,varname,ml_time,0,ml_t3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R T ',trim(ml_fn) |
| varname='U' |
| call getdat(cdfid,varname,ml_time,0,ml_u3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R U ',trim(ml_fn) |
| varname='V' |
| call getdat(cdfid,varname,ml_time,0,ml_v3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R V ',trim(ml_fn) |
| varname='Z' |
| call getdat(cdfid,varname,ml_time,0,ml_z3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R Z ',trim(ml_fn) |
| varname='Q' |
| call getdat(cdfid,varname,ml_time,0,ml_q3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R Q ',trim(ml_fn) |
| varname='PS' |
| call getdat(cdfid,varname,ml_time,0,ml_ps,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R PS ',trim(ml_fn) |
| varname='ORO' |
| call getdat(cdfid,varname,ml_time,0,ml_oro,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ORO ',trim(ml_fn) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Read variables from GEO.MOD |
| call cdfopn(zm_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| varname='T' |
| call getdat(cdfid,varname,zm_time,0,zm_t3,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'R T ',trim(zm_fn) |
| varname='U' |
| call getdat(cdfid,varname,zm_time,0,zm_u3,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'R U ',trim(zm_fn) |
| varname='V' |
| call getdat(cdfid,varname,zm_time,0,zm_v3,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'R V ',trim(zm_fn) |
| varname='P' |
| call getdat(cdfid,varname,zm_time,0,zm_p3,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'R P ',trim(zm_fn) |
| varname='RLAT' |
| call getdat(cdfid,varname,zm_time,0,zm_rlat,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'R RLAT ',trim(zm_fn) |
| varname='RLON' |
| call getdat(cdfid,varname,zm_time,0,zm_rlon,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'R RLON ',trim(zm_fn) |
| |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Read variables from GEO.ORG |
| call cdfopn(or_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| varname='P' |
| call getdat(cdfid,varname,or_time,0,or_p3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R P ',trim(or_fn) |
| varname='U' |
| call getdat(cdfid,varname,or_time,0,or_u3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R U ',trim(or_fn) |
| varname='V' |
| call getdat(cdfid,varname,or_time,0,or_v3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R V ',trim(or_fn) |
| varname='T' |
| call getdat(cdfid,varname,or_time,0,or_t3,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R T ',trim(or_fn) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Initialize the height levels of the Z file |
| do i=1,zm_nx |
| do j=1,zm_ny |
| do k=1,zm_nz |
| zm_z3(i,j,k)=zm_aklay(k) |
| enddo |
| enddo |
| enddo |
| |
| c Calculate 3d pressure field on model levels |
| print*,'C P (ORIGINAL)' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| if (abs(ml_stag(3)+0.5).lt.eps) then |
| ml_p0(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j) |
| else |
| ml_p0(i,j,k)=ml_aklev(k)+ml_bklev(k)*ml_ps(i,j) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c Calculate 3d anomalies due to inversion |
| print*,'C DP (MODIFIED-ORIGINAL)' |
| print*,'C DU (MODIFIED-ORIGINAL)' |
| print*,'C DV (MODIFIED-ORIGINAL)' |
| print*,'C DT (MODIFIED-ORIGINAL)' |
| do k=1,zm_nz |
| do i=1,zm_nx |
| do j=1,zm_ny |
| |
| c P |
| if ( ( abs(zm_p3(i,j,k)-zm_mdv).gt.eps).and. |
| > ( abs(or_p3(i,j,k)-or_mdv).gt.eps) ) |
| > then |
| an_p3(i,j,k) = zm_p3(i,j,k) - or_p3(i,j,k) |
| elseif ( poisson.eq.1 ) then |
| an_p3(i,j,k) = zm_mdv |
| else |
| an_p3(i,j,k) = 0. |
| endif |
| |
| c T |
| if ( ( abs(zm_t3(i,j,k)-zm_mdv).gt.eps).and. |
| > ( abs(or_t3(i,j,k)-or_mdv).gt.eps) ) |
| > then |
| an_t3(i,j,k) = zm_t3(i,j,k) - or_t3(i,j,k) |
| elseif ( poisson.eq.1 ) then |
| an_t3(i,j,k) = zm_mdv |
| else |
| an_t3(i,j,k) = 0. |
| endif |
| |
| c U |
| if ( ( abs(zm_u3(i,j,k)-zm_mdv).gt.eps).and. |
| > ( abs(or_u3(i,j,k)-or_mdv).gt.eps) ) |
| > then |
| an_u3(i,j,k) = zm_u3(i,j,k) - or_u3(i,j,k) |
| elseif ( poisson.eq.1 ) then |
| an_u3(i,j,k) = zm_mdv |
| else |
| an_u3(i,j,k) = 0. |
| endif |
| |
| c V |
| if ( ( abs(zm_v3(i,j,k)-zm_mdv).gt.eps).and. |
| > ( abs(or_v3(i,j,k)-or_mdv).gt.eps) ) |
| > then |
| an_v3(i,j,k) = zm_v3(i,j,k) - or_v3(i,j,k) |
| elseif ( poisson.eq.1 ) then |
| an_v3(i,j,k) = zm_mdv |
| else |
| an_v3(i,j,k) = 0. |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| c ---------------------------------------------------------------- |
| c Get the weight function for the inset (from 1 inside to 0 outside) |
| c ---------------------------------------------------------------- |
| |
| print*,'C BOUNDARY FILTER' |
| |
| c Init the weight function (1 inside, 0 outside) |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| if ( (zm_rlat(i,j).lt. -90.).or. |
| > (zm_rlat(i,j).gt. 90.).or. |
| > (zm_rlon(i,j).lt.-180.).or. |
| > (zm_rlon(i,j).gt. 180.) ) then |
| weight(i,j)=0. |
| else |
| weight(i,j)=1. |
| endif |
| |
| enddo |
| enddo |
| |
| c Get a list of all boundary points |
| nbound=0 |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Get neighbouring points |
| ir=i+1 |
| if (ir.gt.ml_nx) ir=ml_nx |
| il=i-1 |
| if (il.lt. 1) il=1 |
| ju=j+1 |
| if (ju.gt.ml_ny) ju=ml_ny |
| jd=j-1 |
| if (jd.lt. 1) jd=1 |
| |
| c A boundary point has a 0/1 switch |
| if (abs(weight(i,j)-0.).lt.eps) then |
| |
| if ( (abs(weight(ir, j)-1.).lt.eps).or. |
| > (abs(weight(il, j)-1.).lt.eps).or. |
| > (abs(weight(i ,ju)-1.).lt.eps).or. |
| > (abs(weight(i ,jd)-1.).lt.eps).or. |
| > (abs(weight(ir,ju)-1.).lt.eps).or. |
| > (abs(weight(ir,jd)-1.).lt.eps).or. |
| > (abs(weight(il,ju)-1.).lt.eps).or. |
| > (abs(weight(il,jd)-1.).lt.eps).or. |
| > (abs(weight(ir,ju)-1.).lt.eps).or. |
| > (abs(weight(il,ju)-1.).lt.eps).or. |
| > (abs(weight(ir,jd)-1.).lt.eps).or. |
| > (abs(weight(il,jd)-1.).lt.eps) ) then |
| |
| nbound=nbound+1 |
| boundlon(nbound)=ml_xmin+real(i-1)*ml_dx |
| boundlat(nbound)=ml_ymin+real(j-1)*ml_dy |
| |
| endif |
| |
| endif |
| |
| enddo |
| enddo |
| |
| c Get the distance from the subdomain |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| lon=ml_xmin+real(i-1)*ml_dx |
| lat=ml_ymin+real(j-1)*ml_dy |
| |
| distmin=sdis(lon,lat,boundlon(1),boundlat(1)) |
| do k=2,nbound |
| distpos=sdis(lon,lat,boundlon(k),boundlat(k)) |
| if (distpos.lt.distmin) distmin=distpos |
| enddo |
| |
| if ( abs(weight(i,j)-1.).lt.eps) then |
| dist(i,j)=distmin |
| else |
| dist(i,j)=-distmin |
| endif |
| |
| enddo |
| enddo |
| |
| c Set the new weights |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| if (weight(i,j).gt.0.) then |
| weight(i,j)=kink(dist(i,j),transxy) |
| endif |
| |
| enddo |
| enddo |
| |
| c ---------------------------------------------------------------- |
| c Remove MDV regions in input field |
| c ---------------------------------------------------------------- |
| |
| if ( poisson.ne.1 ) goto 120 |
| |
| c Define region for mdv filling |
| do i=1,zm_nx |
| do j=1,zm_ny |
| do k=1,zm_nz |
| |
| c Get neighbour of grid point |
| il = i-1 |
| if (il.lt.1) il=1 |
| ir = i+1 |
| if ( ir.gt.zm_nx ) ir=zm_nx |
| jd = j-1 |
| if (jd.lt.1) jd=1 |
| ju = j+1 |
| if ( ju.gt.zm_ny ) ju=zm_ny |
| |
| c Set flag 2 for boundary and exterior points |
| if ( (abs(zm_rlat(il, j)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(ir, j)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(i , j)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(il,ju)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(ir,ju)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(i ,ju)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(il,jd)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(ir,jd)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(i ,jd)-zm_mdv).lt.eps) ) |
| > then |
| flag_zm(i,j,k) = 2 |
| an_p3(i,j,k) = 0. |
| an_u3(i,j,k) = 0. |
| an_v3(i,j,k) = 0. |
| an_t3(i,j,k) = 0. |
| |
| c Set flag 1 for interior MDV points |
| elseif ( abs(an_p3(i,j,k)-zm_mdv).le.eps ) then |
| flag_zm(i,j,k) = 1 |
| |
| c Set flag 0 for interior valid points |
| else |
| flag_zm(i,j,k) = 0 |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| c Apply mdv filling |
| print*,'C DP (POISSON FILLING)' |
| call mdvfill(an_p3,an_p3,flag_zm,zm_nx,zm_ny,zm_nz,100) |
| print*,'C DT (POISSON FILLING)' |
| call mdvfill(an_t3,an_t3,flag_zm,zm_nx,zm_ny,zm_nz,100) |
| print*,'C DU (POISSON FILLING)' |
| call mdvfill(an_u3,an_u3,flag_zm,zm_nx,zm_ny,zm_nz,100) |
| print*,'C DV (POISSON FILLING)' |
| call mdvfill(an_v3,an_v3,flag_zm,zm_nx,zm_ny,zm_nz,100) |
| |
| c Special treatment: if the number of missing values surpasses 50% |
| c on a level, then the anomaly is imported from the level above |
| do k=zm_nz-1,1,-1 |
| |
| count1 = 0 |
| count0 = 0 |
| do i=1,zm_nx |
| do j= 1,zm_ny |
| if ( flag_zm(i,j,k).eq.1 ) count1 = count1 + 1 |
| if ( flag_zm(i,j,k).eq.0 ) count0 = count0 + 1 |
| enddo |
| enddo |
| |
| if ( count1.gt.count0 ) then |
| print*,'C P (IMPORTING FROM LEVEL ABOVE)',k |
| do i=1,zm_nx |
| do j= 1,zm_ny |
| an_p3(i,j,k) = an_p3(i,j,k+1) |
| flag_zm(i,j,k) = flag_zm(i,j,k+1) |
| enddo |
| enddo |
| print*,'C U (IMPORTING FROM LEVEL ABOVE)',k |
| do i=1,zm_nx |
| do j= 1,zm_ny |
| an_u3(i,j,k) = an_u3(i,j,k+1) |
| flag_zm(i,j,k) = flag_zm(i,j,k+1) |
| enddo |
| enddo |
| |
| print*,'C V (IMPORTING FROM LEVEL ABOVE)',k |
| do i=1,zm_nx |
| do j= 1,zm_ny |
| an_v3(i,j,k) = an_v3(i,j,k+1) |
| flag_zm(i,j,k) = flag_zm(i,j,k+1) |
| enddo |
| enddo |
| |
| print*,'C T (IMPORTING FROM LEVEL ABOVE)',k |
| do i=1,zm_nx |
| do j= 1,zm_ny |
| an_t3(i,j,k) = an_t3(i,j,k+1) |
| flag_zm(i,j,k) = flag_zm(i,j,k+1) |
| enddo |
| enddo |
| |
| endif |
| |
| enddo |
| |
| c Write new fields for tests |
| if ( test.eq.1 ) then |
| |
| call cdfwopn(zm_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 994 |
| |
| isok=0 |
| varname='P_ANO' |
| call check_varok(isok,varname,zm_vnam,zm_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim, |
| > zm_varmin,zm_varmax,zm_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,zm_time,0,an_p3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W P_ANO ',trim(zm_fn) |
| |
| isok=0 |
| varname='T_ANO' |
| call check_varok(isok,varname,zm_vnam,zm_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim, |
| > zm_varmin,zm_varmax,zm_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,zm_time,0,an_t3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W T_ANO ',trim(zm_fn) |
| |
| isok=0 |
| varname='U_ANO' |
| call check_varok(isok,varname,zm_vnam,zm_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim, |
| > zm_varmin,zm_varmax,zm_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,zm_time,0,an_u3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W U_ANO ',trim(zm_fn) |
| |
| isok=0 |
| varname='V_ANO' |
| call check_varok(isok,varname,zm_vnam,zm_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,zm_ndim,zm_mdv,zm_vardim, |
| > zm_varmin,zm_varmax,zm_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,zm_time,0,an_v3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W V_ANO ',trim(zm_fn) |
| |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 994 |
| |
| endif |
| |
| c Exit point for SOR solver |
| 120 continue |
| |
| c ---------------------------------------------------------------- |
| c Change surface pressure and get 3d presure field |
| c ---------------------------------------------------------------- |
| |
| print*,'C PS' |
| |
| c Change surface pressure: based on PV inversion on GEO |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make vertical profile of pressure available |
| n1=0 |
| do k=1,zm_nz |
| n1=n1+1 |
| p1(n1)=an_p3(i,j,k) |
| z1(n1)=zm_z3(i,j,k) |
| enddo |
| |
| if ( n1.ne.0 ) then |
| |
| c Keep the original surface pressure (psfc=0) |
| if ( abs(psfc).lt.eps ) then |
| ml_ps(i,j) = ml_ps(i,j) |
| |
| c Take the full change of surface pressure (psfc=1); |
| c Interpolation/extrapolation with a natural cubic spline |
| elseif ( abs(psfc-1.).lt.eps ) then |
| call spline (z1,p1,n1,1.e30,1.e30,spline_f1) |
| call splint (z1,p1,spline_f1,n1,ml_oro(i,j),hh) |
| ml_ps(i,j)=ml_ps(i,j) + hh*weight(i,j) |
| |
| c Only take a fractional change of surface pressure |
| c Interpolation/extrapolation with a natural cubic spline |
| else |
| call spline (z1,p1,n1,1.e30,1.e30,spline_f1) |
| call splint (z1,p1,spline_f1,n1,ml_oro(i,j),hh) |
| ml_ps(i,j)=ml_ps(i,j) + hh*psfc*weight(i,j) |
| |
| endif |
| |
| endif |
| |
| enddo |
| enddo |
| |
| c Calculate 3d pressure field on model levels |
| print*,'C P' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| if (abs(ml_stag(3)+0.5).lt.eps) then |
| ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j) |
| else |
| ml_p3(i,j,k)=ml_aklev(k)+ml_bklev(k)*ml_ps(i,j) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c ---------------------------------------------------------------- |
| c Get T,U,V at the new pressure levels, based on the original P file |
| c ---------------------------------------------------------------- |
| |
| c Interpolate T from original P file to new pressure levels |
| print*,'C T (ORIGINAL)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| n1=0 |
| do k=ml_nz,1,-1 |
| n1=n1+1 |
| f1(n1)=ml_t3(i,j,k) |
| p1(n1)=ml_p0(i,j,k) |
| enddo |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_t3(i,j,k)=hh |
| enddo |
| |
| enddo |
| enddo |
| |
| c Interpolate U from original P file to new pressure levels |
| print*,'C U (ORIGINAL)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| n1=0 |
| do k=ml_nz,1,-1 |
| n1=n1+1 |
| f1(n1)=ml_u3(i,j,k) |
| p1(n1)=ml_p0(i,j,k) |
| enddo |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_u3(i,j,k)=hh |
| enddo |
| |
| enddo |
| enddo |
| |
| c Interpolate V from original P file to new pressure levels |
| print*,'C V (ORIGINAL)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| n1=0 |
| do k=ml_nz,1,-1 |
| n1=n1+1 |
| f1(n1)=ml_v3(i,j,k) |
| p1(n1)=ml_p0(i,j,k) |
| enddo |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_v3(i,j,k)=hh |
| enddo |
| |
| enddo |
| enddo |
| |
| c ---------------------------------------------------------------- |
| c Add T,U,V anomalies at the new pressure levels |
| c ---------------------------------------------------------------- |
| |
| c Change temperature field |
| print*,'C T (ANOMALY)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make vertical profile of temperature available |
| n1=0 |
| pmax=-100. |
| pmin=2000. |
| do k=zm_nz,1,-1 |
| if ((abs(an_t3(i,j,k)-zm_mdv).gt.eps).and. |
| > (zm_z3(i,j,k).gt.ml_oro(i,j)).and. |
| > (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then |
| n1=n1+1 |
| f1(n1)=an_t3(i,j,k) |
| p1(n1)=zm_p3(i,j,k) |
| if (p1(n1).lt.pmin) pmin=p1(n1) |
| if (p1(n1).gt.pmax) pmax=p1(n1) |
| endif |
| enddo |
| pmin=pmin-dpextra*pmin |
| pmax=pmax+dpextra+pmax |
| |
| c Cubic spline interpolation |
| if (n1.gt.0) then |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| if ( (ml_p3(i,j,k).gt.pmin).and. |
| > (ml_p3(i,j,k).lt.pmax) ) then |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_t3(i,j,k)=ml_t3(i,j,k) + hh*weight(i,j) |
| endif |
| enddo |
| endif |
| |
| enddo |
| enddo |
| |
| c Change zonal velocity field |
| print*,'C U' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make vertical profile of zonal velocity available |
| n1=0 |
| pmax=-100. |
| pmin=2000. |
| do k=zm_nz,1,-1 |
| if ((abs(an_u3(i,j,k)-zm_mdv).gt.eps).and. |
| > (zm_z3(i,j,k).gt.ml_oro(i,j)).and. |
| > (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then |
| n1=n1+1 |
| f1(n1)=an_u3(i,j,k) |
| p1(n1)=zm_p3(i,j,k) |
| if (p1(n1).lt.pmin) pmin=p1(n1) |
| if (p1(n1).gt.pmax) pmax=p1(n1) |
| endif |
| enddo |
| pmin=pmin-dpextra*pmin |
| pmax=pmax+dpextra*pmax |
| |
| c Cubic spline interpolation |
| if (n1.gt.0) then |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| if ( (ml_p3(i,j,k).gt.pmin).and. |
| > (ml_p3(i,j,k).lt.pmax) ) then |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_u3(i,j,k)=ml_u3(i,j,k) + hh*weight(i,j) |
| endif |
| enddo |
| endif |
| |
| enddo |
| enddo |
| |
| c Change meridional velocity field |
| print*,'C V' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make vertical profile of zonal velocity available |
| n1=0 |
| pmax=-100. |
| pmin=2000. |
| do k=zm_nz,1,-1 |
| if ((abs(an_v3(i,j,k)-zm_mdv).gt.eps).and. |
| > (zm_z3(i,j,k).gt.ml_oro(i,j)).and. |
| > (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then |
| n1=n1+1 |
| f1(n1)=an_v3(i,j,k) |
| p1(n1)=zm_p3(i,j,k) |
| if (p1(n1).lt.pmin) pmin=p1(n1) |
| if (p1(n1).gt.pmax) pmax=p1(n1) |
| endif |
| enddo |
| pmin=pmin-dpextra*pmin |
| pmax=pmax+dpextra*pmax |
| |
| c Cubic spline interpolation |
| if (n1.gt.0) then |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| if ( (ml_p3(i,j,k).gt.pmin).and. |
| > (ml_p3(i,j,k).lt.pmax) ) then |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_v3(i,j,k)=ml_v3(i,j,k) + hh*weight(i,j) |
| endif |
| enddo |
| endif |
| |
| enddo |
| enddo |
| |
| c --------------------------------------------------------------------- |
| c Change geopotential height |
| c --------------------------------------------------------------------- |
| |
| c Interpolate Z from original P file to new pressure levels |
| print*,'C Z (ORIGINAL)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| n1=0 |
| do k=ml_nz,1,-1 |
| n1=n1+1 |
| f1(n1)=ml_z3(i,j,k) |
| p1(n1)=ml_p0(i,j,k) |
| enddo |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| ml_z3(i,j,k)=hh |
| enddo |
| |
| enddo |
| enddo |
| |
| c Calculate 3d virtual temperature |
| print*,'C TV' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)* |
| > (1.+kappa*ml_q3(i,j,k)) |
| enddo |
| enddo |
| enddo |
| |
| c Add geopotential height anomaly: first, the pressure anomaly is |
| c interpolated to the new position, then it is transformed into |
| c a correction of geopotential height with the hydrostatic equation |
| print*,'C DZ (MODIFIED -ORIGINAL)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make vertical profile of pressure available |
| n1=0 |
| pmax=-100. |
| pmin=2000. |
| do k=zm_nz,1,-1 |
| if ((abs(an_p3(i,j,k)-zm_mdv).gt.eps).and. |
| > (zm_z3(i,j,k).gt.ml_oro(i,j)).and. |
| > (abs(zm_p3(i,j,k)-zm_mdv).gt.eps)) then |
| n1=n1+1 |
| f1(n1)=an_p3(i,j,k) |
| p1(n1)=zm_p3(i,j,k) |
| if (p1(n1).lt.pmin) pmin=p1(n1) |
| if (p1(n1).gt.pmax) pmax=p1(n1) |
| endif |
| enddo |
| pmin=pmin-dpextra*pmin |
| pmax=pmax+dpextra*pmax |
| |
| c Cubic spline interpolation and conversion dp -> dz |
| if (n1.gt.0) then |
| call spline (p1,f1,n1,1.e30,1.e30,spline_f1) |
| do k=1,ml_nz |
| if ( (ml_p3(i,j,k).gt.pmin).and. |
| > (ml_p3(i,j,k).lt.pmax) ) then |
| call splint (p1,f1,spline_f1,n1,ml_p3(i,j,k),hh) |
| hh = -rdg * ml_tv3(i,j,k) * hh/ml_p3(i,j,k) |
| ml_z3(i,j,k)=ml_z3(i,j,k) + hh*weight(i,j) |
| endif |
| enddo |
| endif |
| |
| enddo |
| enddo |
| |
| c ---------------------------------------------------------------- |
| c Remove unrealistic values from final fields |
| c ---------------------------------------------------------------- |
| |
| if ( poisson.ne.1 ) goto 120 |
| |
| c Define region for mdv filling |
| count1 = 0 |
| |
| do i=1,ml_nx |
| do j=1,ml_ny |
| do k=1,ml_nz |
| |
| c Get neighbour of grid point |
| il = i-1 |
| if (il.lt.1) il=1 |
| ir = i+1 |
| if ( ir.gt.ml_nx ) ir=ml_nx |
| jd = j-1 |
| if (jd.lt.1) jd=1 |
| ju = j+1 |
| if ( ju.gt.ml_ny ) ju=ml_ny |
| |
| c Set flag 2 for boundary and exterior points |
| if ( (abs(zm_rlat(il, j)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(ir, j)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(i , j)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(il,ju)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(ir,ju)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(i ,ju)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(il,jd)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(ir,jd)-zm_mdv).lt.eps).or. |
| > (abs(zm_rlat(i ,jd)-zm_mdv).lt.eps) ) |
| > then |
| flag_ml(i,j,k) = 2 |
| |
| c Set flag 1 for interior unphysical points |
| elseif ( ( abs(ml_t3(i,j,k)).gt.500. ).or. |
| > ( abs(ml_u3(i,j,k)).gt.500. ).or. |
| > ( abs(ml_v3(i,j,k)).gt.500. ) ) |
| > then |
| flag_ml(i,j,k) = 1 |
| count1 = count1 + 1 |
| |
| c Set flag 0 for interior valid points |
| else |
| flag_ml(i,j,k) = 0 |
| endif |
| |
| enddo |
| enddo |
| enddo |
| |
| print*,'C MASK UNREALISTIC VALUES FOR T,U,V',count1 |
| |
| c Apply mdv filling |
| print*,'C T (POISSON FILLING)' |
| call mdvfill(ml_t3,ml_t3,flag_ml,ml_nx,ml_ny,ml_nz,10) |
| print*,'C U (POISSON FILLING)' |
| call mdvfill(ml_u3,ml_u3,flag_ml,ml_nx,ml_ny,ml_nz,10) |
| print*,'C V (POISSON FILLING)' |
| call mdvfill(ml_v3,ml_v3,flag_ml,ml_nx,ml_ny,ml_nz,10) |
| |
| c ---------------------------------------------------------------- |
| c Write modified fields to P file |
| c ---------------------------------------------------------------- |
| |
| c Open output file |
| call cdfwopn(ml_fn,cdfid,ierr) |
| if (ierr.ne.0) goto 994 |
| |
| c Write surface pressure |
| isok=0 |
| varname='PS' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| ml_vardim(3)=3 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| ml_vardim(3)=4 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_ps,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W PS ',trim(ml_fn) |
| |
| c Write temperature |
| isok=0 |
| varname='T' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_t3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W T ',trim(ml_fn) |
| |
| c Write zonal wind |
| isok=0 |
| varname='U' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_u3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W U ',trim(ml_fn) |
| |
| c Write meridional wind |
| isok=0 |
| varname='V' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_v3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W V ',trim(ml_fn) |
| |
| c Write geopotential height |
| isok=0 |
| varname='Z' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_z3,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W Z ',trim(ml_fn) |
| |
| c Filter matrix (only in test mode) |
| if (test.eq.1) then |
| isok=0 |
| varname='DIST' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| ml_vardim(3)=1 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,ml_time,1,dist,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W DIST ',trim(ml_fn) |
| |
| isok=0 |
| varname='WEIGHT' |
| call check_varok(isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| ml_vardim(3)=1 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 994 |
| endif |
| call putdat(cdfid,varname,ml_time,1,weight,ierr) |
| if (ierr.ne.0) goto 994 |
| print*,'W WEIGHT ',trim(ml_fn) |
| |
| endif |
| |
| c Close output file |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 994 |
| |
| c ---------------------------------------------------------------- |
| c Exception handling |
| c ---------------------------------------------------------------- |
| |
| stop |
| |
| 998 print*,'Problem with input netcdf file ',trim(ml_fn) |
| stop |
| |
| 997 print*,'Problem with input netcdf file ',trim(zm_fn) |
| stop |
| |
| 994 print*,'Problem with output netcdf file ',trim(ml_fn) |
| 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 Natural cubic spline /from Numerical Recipes) |
| c ------------------------------------------------------------- |
| |
| SUBROUTINE spline(x,y,n,yp1,ypn,y2) |
| INTEGER n,NMAX |
| REAL yp1,ypn,x(n),y(n),y2(n) |
| PARAMETER (NMAX=500) |
| INTEGER i,k |
| REAL p,qn,sig,un,u(NMAX) |
| if (yp1.gt..99e30) then |
| y2(1)=0. |
| u(1)=0. |
| else |
| y2(1)=-0.5 |
| u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) |
| endif |
| do 11 i=2,n-1 |
| sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) |
| p=sig*y2(i-1)+2. |
| y2(i)=(sig-1.)/p |
| u(i)=(6.*((y(i+1)-y(i))/(x(i+ |
| *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig* |
| *u(i-1))/p |
| 11 continue |
| if (ypn.gt..99e30) then |
| qn=0. |
| un=0. |
| else |
| qn=0.5 |
| un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) |
| endif |
| y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) |
| do 12 k=n-1,1,-1 |
| y2(k)=y2(k)*y2(k+1)+u(k) |
| 12 continue |
| return |
| END |
| |
| |
| SUBROUTINE splint(xa,ya,y2a,n,x,y) |
| INTEGER n |
| REAL x,y,xa(n),y2a(n),ya(n) |
| INTEGER k,khi,klo |
| REAL a,b,h |
| klo=1 |
| khi=n |
| 1 if (khi-klo.gt.1) then |
| k=(khi+klo)/2 |
| if(xa(k).gt.x)then |
| khi=k |
| else |
| klo=k |
| endif |
| goto 1 |
| endif |
| h=xa(khi)-xa(klo) |
| if (h.eq.0.) pause 'bad xa input in splint' |
| a=(xa(khi)-x)/h |
| b=(x-xa(klo))/h |
| y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** |
| *2)/6. |
| return |
| END |
| |
| |
| c --------------------------------------------------------- |
| c Spherical distance between two lat/lon points |
| c --------------------------------------------------------- |
| |
| real function sdis(xp,yp,xq,yq) |
| c |
| c calculates spherical distance (in km) between two points given |
| c by their spherical coordinates (xp,yp) and (xq,yq), respectively. |
| c |
| real re |
| parameter (re=6370.) |
| real xp,yp,xq,yq,arg |
| |
| arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq) |
| if (arg.lt.-1.) arg=-1. |
| if (arg.gt.1.) arg=1. |
| sdis=re*acos(arg) |
| |
| end |
| |
| C ----------------------------------------------------------------- |
| c Rene's KINK function (for smoothing at bounadries) |
| C ----------------------------------------------------------------- |
| |
| real function kink (x,a) |
| |
| implicit none |
| |
| c declaration of parameters |
| real x,a |
| |
| c parameters |
| real pi |
| parameter (pi=3.1415926535) |
| |
| if (x.lt.0.) then |
| kink=0. |
| elseif (x.gt.a) then |
| kink=1. |
| else |
| kink=(1.+cos(pi*(x-a)/a))/2. |
| endif |
| |
| return |
| end |
| |
| C ----------------------------------------------------------------- |
| c Rene's KINK function (for smoothing at bounadries) |
| C ----------------------------------------------------------------- |
| |
| subroutine mdvfill (out,inp,flag,nx,ny,nz,maxiter) |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real inp (nx,ny,nz) |
| real out (nx,ny,nz) |
| integer flag(nx,ny,nz) |
| integer maxiter |
| |
| c Parameters |
| real omega ! Omega fopr SOR |
| parameter (omega=1.5) |
| |
| c Auxiliary variables |
| integer i,j,k |
| integer iter |
| real tmp0(nx,ny,nz) |
| real tmp1(nx,ny,nz) |
| integer il,ir,ju,jd |
| integer count |
| real mean(nz) |
| |
| c Calculate mean of variable for all levels |
| do k=1,nz |
| mean(k) = 0. |
| count = 0 |
| do i=1,nx |
| do j=1,ny |
| if ( flag(i,j,k).eq.0 ) then |
| count = count + 1 |
| mean(k) = mean(k) + inp(i,j,k) |
| endif |
| enddo |
| enddo |
| if ( count.ne.0 ) then |
| mean(k) = mean(k)/real(count) |
| else |
| mean(k) = 0. |
| endif |
| enddo |
| |
| c Create first guess |
| do k=1,nz |
| do i=1,nx |
| do j=1,ny |
| if ( flag(i,j,k).eq.0 ) then |
| tmp0(i,j,k) = inp(i,j,k) |
| else |
| tmp0(i,j,k) = mean(k) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c SOR iterations |
| iter = 0 |
| 122 continue |
| |
| c Loop over all points |
| do k=1,nz |
| do i=1,nx |
| do j=1,ny |
| |
| c Apply the updating only for specified points |
| if ( flag(i,j,k).ne.1 ) goto 121 |
| |
| c Get neighbouring points - no handling of periodic domains! |
| il = i-1 |
| if (il.lt.1) il=1 |
| ir = i+1 |
| if ( ir.gt.nx ) ir=nx |
| jd = j-1 |
| if (jd.lt.1) jd=1 |
| ju = j+1 |
| if ( ju.gt.ny ) ju=ny |
| |
| c Update field |
| tmp1(i,j,k) = 0.25 * ( tmp0(il,j,k) + tmp0(ir,j,k) + |
| > tmp0(i,ju,k) + tmp0(i,jd,k) ) |
| |
| tmp0(i,j,k) = omega * tmp1(i,j,k) + |
| > (1. - omega) * tmp1(i,j,k) |
| |
| c Exit point for loop |
| 121 continue |
| |
| enddo |
| enddo |
| enddo |
| |
| c Decide whether further iterations are needed |
| iter = iter + 1 |
| if ( iter.lt.maxiter) goto 122 |
| |
| c Save output |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ( flag(i,j,k).eq.1 ) then |
| out(i,j,k) = tmp0(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| end |
| |
| c |
| Property changes: |
|---|
| Added: svn:executable |
|---|