program p2cinc use mkinds use mncio use mdates use mabio use mutil implicit none integer(i4),parameter :: src_nx=525,src_ny=385,src_nz=41,src_nt=1 integer(i4),parameter :: dst_nx=2101,dst_ny=1537,dst_nz=41,dst_nt=1 character(len=cl) :: srcFile,dstFile real(r4),allocatable :: src_fld(:,:),dst_fld(:,:),depth(:,:) real(r4),allocatable :: src_lon(:,:),src_lat(:,:),dst_lon(:,:),dst_lat(:,:) real(r4) :: rlmask(src_nx,src_ny) logical(bl) :: lmask(src_nx,src_ny) integer(i4) :: iw(src_nx,src_ny) real(r4),allocatable :: ftem(:,:,:),fsal(:,:,:),fuvl(:,:,:),fvvl(:,:,:),fthk(:,:,:) real(r4),allocatable :: fuba(:,:),fvba(:,:),fssh(:,:) real(r4),allocatable :: inc_tem(:,:,:),inc_sal(:,:,:),inc_uvl(:,:,:),inc_vvl(:,:,:),inc_thk(:,:,:) logical,parameter :: ltides=.false. !!Alex ! Internal arrays INTEGER(i4) :: xmap(dst_nx,dst_ny,1,4) INTEGER(i4) :: ymap(dst_nx,dst_ny,1,4) INTEGER(i4) :: mask(dst_nx,dst_ny) REAL(r4) :: wgts(dst_nx,dst_ny,1,4) real(r4),allocatable :: xax(:),yax(:),map(:,:,:,:),mlayer(:) ! internal variables integer(i4) :: den,i,j,k,l,fid,offset,iif,iil,jjf,jjl,nf,iwp,npass,sver character(:),allocatable :: fname,filename,in_rst_file,in_inc_file,out_inc_file,outFile,mapfile integer(i4) :: vtype,iunit_in,iunit_out integer(i4) :: vdims(4),vdim(1),vdims3(3),vdims2(2) INTEGER(i4) :: ii INTEGER(i4) :: jj REAL(r4) :: ptlon REAL(r4) :: ptlat REAL(r4) :: grdi REAL(r4) :: grdj REAL(r4) :: r REAL(r4) :: s REAL(r4) :: xsum,grdt,val,thbase,coord REAL(r4),PARAMETER :: spval=2.00**100 character(10) :: cdtg real(r4) :: tstep,mday,sday,hmin,hmax CHARACTER(len=2) :: fn integer(i4) :: iyear,iyday,ihour,nleap character(4) :: cyear character(3) :: cday character(3) :: chour character(4) :: clevel character(22) :: sflnm character(:),allocatable :: gomdataloc,gfsdataloc,mwvdataloc real(r4) :: time,atime(1) character(cl) :: buffer CHARACTER(len=120), DIMENSION(5) :: preambl CHARACTER(len=120) :: cline character(len=120) :: ctitle(4) integer(i4) :: yrflag,iversn,sigver,iexpt,nstep,iunit integer(i4) :: iv1 integer(i4) :: iv2 real(r4) :: rv3 real(r4) :: rv4 character(len=8) :: vname print*,'... Entering p2cinc ....' !bfile headers ctitle(1)="TSIS Analysis for incremental update" ctitle(2)="GDEM3 Jan. init; KPP; SeaWiFS KPAR; e-loan ice; A=20; Smag=.05; veldf4=.02;" ctitle(3)="S-Z(14-18):dp00/f/x/i=3m/1.18/120m/1m;ds=0.5m/1.18/75m; var den: Med+Black Seas;" ctitle(4)="Dummy Info in the above two lines" ! iversn=22 iexpt=191 yrflag=3 nstep=2999100 thbase=34.0 iunit_in=11 iunit_out=12 call getarg(1,buffer) in_rst_file=trim(buffer) call getarg(2,buffer) in_inc_file=trim(buffer) call getarg(3,buffer) out_inc_File=trim(buffer) call getarg(4,cdtg) call CDATE2WNDAY(time,cdtg) !time=time-0.25 print*,'... Arguments p2cinc ....' allocate(src_fld(src_nx,src_ny)) allocate(dst_fld(dst_nx,dst_ny)) allocate(depth(dst_nx,dst_ny)) allocate(src_lon(src_nx,src_ny)) allocate(src_lat(src_nx,src_ny)) allocate(dst_lon(dst_nx,dst_ny)) allocate(dst_lat(dst_nx,dst_ny)) allocate(ftem(dst_nx,dst_ny,dst_nz)) allocate(fsal(dst_nx,dst_ny,dst_nz)) allocate(fthk(dst_nx,dst_ny,dst_nz)) allocate(fuvl(dst_nx,dst_ny,dst_nz)) allocate(fvvl(dst_nx,dst_ny,dst_nz)) allocate(fuba(dst_nx,dst_ny)) allocate(fvba(dst_nx,dst_ny)) allocate(fssh(dst_nx,dst_ny)) allocate(inc_uvl(src_nx,src_ny,src_nz)) allocate(inc_vvl(src_nx,src_ny,src_nz)) allocate(inc_thk(src_nx,src_ny,src_nz)) allocate(inc_tem(src_nx,src_ny,src_nz)) allocate(inc_sal(src_nx,src_ny,src_nz)) allocate(mlayer(dst_nz)) fname="p_gridinfo.nc" CALL nciopn(fname,fid) CALL nciorv(fname,fid,"mplon",src_lon) CALL nciorv(fname,fid,"mplat",src_lat) CALL nciorv(fname,fid,"target_densities",mlayer) CALL nciocl(fname,fid) fname="c_gridinfo.nc" CALL nciopn(fname,fid) CALL nciorv(fname,fid,"mplon",dst_lon) CALL nciorv(fname,fid,"mplat",dst_lat) CALL nciorv(fname,fid,"mdepth",depth) CALL nciocl(fname,fid) where(dst_lon>180.0) dst_lon=dst_lon-360.0 k=1 DO jj=1,dst_ny DO ii=1,dst_nx ptlon=dst_lon(ii,jj) ptlat=dst_lat(ii,jj) CALL ll2ij(src_lon(:,1),src_lat(1,:),ptlon,ptlat,src_nx,src_ny,grdi,grdj) CALL setwgts(src_nx,src_ny,grdi,grdj,xmap(ii,jj,k,:),ymap(ii,jj,k,:),wgts(ii,jj,k,:)) IF(sum(wgts(ii,jj,k,:)) > 1.01) THEN WRITE(*,*) "Weights Error: ", ii,jj STOP ENDIF xsum=xsum + sum(wgts(ii,jj,k,:)) ENDDO ENDDO WRITE(*,*) "Min values of mapping coodinates in src grid: ", src_lat(1,MINVAL(ymap(:,:,k,:))),src_lon(MINVAL(xmap(:,:,k,:)),1) WRITE(*,*) "Max values of mapping coodinates in src grid: ", src_lat(1,MAXVAL(ymap(:,:,k,:))),src_lon(MAXVAL(xmap(:,:,k,:)),1) ! load restart iunit=10 call zaiost(dst_nx,dst_ny) call zaiopn(in_rst_file,iunit) ! skip nz fields and read nz uvl(xward) fields at the second time level DO k=1,src_nz CALL zaiosk(iunit) ENDDO DO k=1,src_nz CALL zaiord(iunit,fuvl(:,:,k),hmin,hmax) ENDDO ! skip nz fields and read nz vvl(yward) fields at the second time level DO k=1,src_nz CALL zaiosk(iunit) ENDDO DO k=1,src_nz CALL zaiord(iunit,fvvl(:,:,k),hmin,hmax) ENDDO ! skip nz fields and read nz thickness fields at the second time level DO k=1,src_nz CALL zaiosk(iunit) ENDDO DO k=1,src_nz CALL zaiord(iunit,fthk(:,:,k),hmin,hmax) ENDDO ! skip nz fields and read nz tem fields at the second time level DO k=1,src_nz CALL zaiosk(iunit) ENDDO DO k=1,src_nz CALL zaiord(iunit,ftem(:,:,k),hmin,hmax) ENDDO ! skip nz fields and read nz sal fields at the second time level DO k=1,src_nz CALL zaiosk(iunit) ENDDO DO k=1,src_nz CALL zaiord(iunit,fsal(:,:,k),hmin,hmax) ENDDO !! skip tidal fields !!Alex if (ltides) then DO k=1,2*49 CALL zaiosk(iunit) ENDDO endif ! now read 2d fields; get the mid-time level ! ubaro CALL zaiosk(iunit) CALL zaiord(iunit,fuba(:,:),hmin,hmax) ! mid-time ubaro CALL zaiosk(iunit) ! vbaro CALL zaiosk(iunit) CALL zaiord(iunit,fvba(:,:),hmin,hmax) ! mid-time vbaro CALL zaiosk(iunit) ! Close .ab files call zaiocl(iunit) ! read inc archv fileName=in_inc_file iunit=11 ! Initialize hycom style IO CALL zaiost(src_nx,src_ny) CALL zaiopn(FileName,iunit) ! First read 2d fields skipping over some fields CALL zaiosk(iunit) ! skip montg !CALL zaiosk(iunit) ! skip ssh CALL zaiord(iunit,fssh(:,:),hmin,hmax) ! read ssh !if(lsteric) then !CALL zaiosk(iunit) ! skip steric !endif CALL zaiosk(iunit) ! skip surflx CALL zaiosk(iunit) ! skip salflx CALL zaiosk(iunit) ! skip bl_dpth CALL zaiosk(iunit) ! skip mx_dpth CALL zaiosk(iunit) ! skip ubaro CALL zaiosk(iunit) ! skip vbaro !CALL zaiord(iunit,fmld(:,:),hmin,hmax) ! read mx_dpth !CALL zaiord(iunit,fuba(:,:),hmin,hmax) ! read ubaro !CALL zaiord(iunit,fvba(:,:),hmin,hmax) ! read vbaro ! Next read layered fields do k=1,src_nz CALL zaiord(iunit,inc_uvl(:,:,k),hmin,hmax) ! read xward velocities CALL zaiord(iunit,inc_vvl(:,:,k),hmin,hmax) ! read yward velocities CALL zaiord(iunit,inc_thk(:,:,k),hmin,hmax) ! read thickness CALL zaiord(iunit,inc_tem(:,:,k),hmin,hmax) ! read pot. temperature CALL zaiord(iunit,inc_sal(:,:,k),hmin,hmax) ! read salinity enddo ! Close .ab files call zaiocl(iunit) do k=1,1 !!Alex remove comment on remapping uvl and vvl !!Alex add the masking of high values of uvl and vvl !uvl ! dst_fld=0. ! call remap(dst_fld,inc_uvl(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) ! fuvl(:,:,k)=fuvl(:,:,k)+ dst_fld ! where(abs(fuvl(:,:,k))>5.0) fuvl(:,:,k)=0.0 !!vvl ! dst_fld=0. ! call remap(dst_fld,inc_vvl(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) ! fvvl(:,:,k)=fvvl(:,:,k)+ dst_fld ! where(abs(fvvl(:,:,k))>5.0) fvvl(:,:,k)=0.0 !tem dst_fld=0. call remap(dst_fld,inc_tem(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) ftem(:,:,k)=ftem(:,:,k)+ dst_fld ! sal dst_fld=0. call remap(dst_fld,inc_sal(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) fsal(:,:,k)=fsal(:,:,k)+ dst_fld enddo !thk do k=1,dst_nz !!Alex remove comment on remapping uvl and vvl !!Alex add the masking of high values of uvl and vvl !uvl dst_fld=0. call remap(dst_fld,inc_uvl(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) fuvl(:,:,k)=fuvl(:,:,k)+ dst_fld where(abs(fuvl(:,:,k))>5.0) fuvl(:,:,k)=0.0 !vvl dst_fld=0. call remap(dst_fld,inc_vvl(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) fvvl(:,:,k)=fvvl(:,:,k)+ dst_fld where(abs(fvvl(:,:,k))>5.0) fvvl(:,:,k)=0.0 dst_fld=0. call remap(dst_fld,inc_thk(:,:,k),xmap(:,:,1,:),ymap(:,:,1,:),wgts(:,:,1,:),src_nx,src_ny,dst_nx,dst_ny) ! print *, fthk(340,140,k),fthk(340,140,k)+ dst_fld(340,140) fthk(:,:,k)=fthk(:,:,k)+ dst_fld ! print *, minval(dst_fld)/9086.,maxval(dst_fld)/9806.,minval(fthk(:,:,k))/9806.,maxval(fthk(:,:,k),mask=fthk(:,:,k)<2.00**99)/9806. enddo ! clean dp do j=1,dst_ny do i=1,dst_nx if(depth(i,j)<2.00**99) then ! clean up any negative thickness 1st pass do k = 1, dst_nz-1 fthk(i,j,k+1) =fthk(i,j,k+1)+min(0.0,fthk(i,j,k)) fthk(i,j,k) =max(fthk(i,j,k),0.0) end do ! clean up any negative thicknes 2nd pass do k = dst_nz, 2, -1 fthk(i,j,k-1) = fthk(i,j,k-1)+min(0.0,fthk(i,j,k)) fthk(i,j,k) = max(fthk(i,j,k),0.0) end do ! end post procesing else fthk(i,j,:)=0.0 endif enddo enddo !do k=1,dst_nz ! print *, minval(fthk(:,:,k))/9806.,maxval(fthk(:,:,k),mask=fthk(:,:,k)<2.00**99)/9806. !enddo fileName=out_inc_file iunit=11 sver=6 ! write out the data call zaiost(dst_nx,dst_ny) call zaiopn(fileName,iunit) write(iunit,116) ctitle,iversn,iexpt,yrflag,dst_nx,dst_ny call zaiowr(iunit,fssh,hmin,hmax) WRITE (iunit,117) 'montg1 ',nstep,time,sver,thbase,hmin,hmax call zaiowr(iunit,fssh,hmin,hmax) WRITE (iunit,117) 'srfhgt ',nstep,time,0,0,hmin,hmax !call zaiowr(iunit,dummy,hmin,hmax) !WRITE (iunit,117) 'steric ',nstep,time,0,0,hmin,hmax call zaiowr(iunit,fssh,hmin,hmax) WRITE (iunit,117) 'surflx ',nstep,time,0,0,hmin,hmax call zaiowr(iunit,fssh,hmin,hmax) WRITE (iunit,117) 'salflx ',nstep,time,0,0,hmin,hmax call zaiowr(iunit,fssh,hmin,hmax) WRITE (iunit,117) 'bl_dpth ',nstep,time,0,0,hmin,hmax call zaiowr(iunit,fssh,hmin,hmax) WRITE (iunit,117) 'mix_dpth ',nstep,time,0,0,hmin,hmax call zaiowr(iunit,fuba,hmin,hmax) WRITE (iunit,117) 'u_btrop ',nstep,time,0,0,hmin,hmax call zaiowr(iunit,fvba,hmin,hmax) WRITE (iunit,117) 'v_btrop ',nstep,time,0,0,hmin,hmax ! layer loop DO k=1,dst_nz coord=mlayer(k) CALL zaiowr(iunit,fuvl(:,:,k),hmin,hmax) print *,"fuvl: ", k, minval(fuvl(:,:,k)),maxval(fuvl(:,:,k)) WRITE (iunit,117) 'u-vel. ',nstep,time,k,coord,hmin,hmax CALL zaiowr(iunit,fvvl(:,:,k),hmin,hmax) print *,"fvvl: ", k, minval(fvvl(:,:,k)),maxval(fvvl(:,:,k)) WRITE (iunit,117) 'v-vel. ',nstep,time,k,coord,hmin,hmax CALL zaiowr(iunit,fthk(:,:,k),hmin,hmax) WRITE (iunit,117) 'thknss ',nstep,time,k,coord,hmin,hmax CALL zaiowr(iunit,ftem(:,:,k),hmin,hmax) WRITE (iunit,117) 'temp ',nstep,time,k,coord,hmin,hmax CALL zaiowr(iunit,fsal(:,:,k),hmin,hmax) WRITE (iunit,117) 'salin ',nstep,time,k,coord,hmin,hmax ENDDO call zaiocl(iunit) 116 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & 'field time step model day',' k dens min max') 117 FORMAT (a8,' =',i11,f11.3,i3,f6.2,1p2e16.7) end program p2cinc