module margo use mkinds use mdates use mncio use mutil use mseawater use mgdem use mgrid use mhycom use miotsis use mpobs implicit none contains subroutine argor2tsidb(cdtg) implicit none character(:),allocatable :: argoFile,dbloc,atlFile,pacFile,indFile,rfile character(10) :: cdtg real(r4) :: wday,wsday real(r4) :: lonmn,lonmx,latmn,latmx,plon(1),plat(1) integer :: ncid,nprof,nlvl,vid,fid real(r4),allocatable :: lon(:),lat(:) integer(i4),parameter :: ntzlevs=52 real(r4) :: tzlevels(ntzlevs) integer :: n_outliers,ninv,i,j,ncount,k,kk,stat,ngdem,nvalid,vcount,fcount,ii,jj,ndups,nwithin,nuniq real(r4),allocatable :: pres(:,:),temp(:,:),saln(:,:) character(:),allocatable :: pres_qc,temp_qc,saln_qc,juld_qc,pos_qc real(r4),allocatable :: pres_qc_wod(:),temp_qc_wod(:),saln_qc_wod(:) integer(i4),allocatable :: valid(:) real(r4),allocatable :: pres1d(:),temp1d(:),saln1d(:),pres1d_qc(:),temp1d_qc(:),saln1d_qc(:) real(r4),allocatable :: ftemp1d(:),fsaln1d(:),ptemp(:) character(4) :: cnprof character(:),allocatable :: FileName character(:), allocatable :: vname integer(i4) :: vtype integer(i4) :: vdims3(3),vdims2(2),vdim(1) integer(i4),allocatable :: xax(:),zax(:) real(r8),allocatable :: juldr8(:) real(r4) :: date_offset logical(bl) :: lexist,lrexist,laexist,lpexist,liexist,dt_argo,rt_argo integer,parameter :: nprofmx=10000 real(r4),allocatable :: nlon(:),nlat(:),nptemp(:,:),nsaln(:,:),ptime(:) real(r4),allocatable :: slon(:),slat(:) integer(i4) :: snx,sny,kmn,kmx integer(i4),allocatable :: indx(:,:),akmx(:,:),prof_kmx(:),prof_valid(:),prof_num(:) real(r4),allocatable :: vtemp(:,:),vsaln(:,:),val(:),vtime(:,:),err(:),verr(:,:) real(r4),allocatable :: flon(:),flat(:),ftime(:,:),ftemp(:,:,:),fsaln(:,:,:),fdens(:,:,:),fuvel(:,:,:),fvvel(:,:,:),prof_dens(:,:) real(r4),allocatable :: tclim(:),sclim(:),tstd(:),sstd(:) real(r4),allocatable :: atclim(:,:),asclim(:,:),atstd(:,:),asstd(:,:) !real(r4),allocatable :: pl(:),rl(:),sl(:),tl(:),ul(:),vl(:),dp(:),terr(:),serr(:),derr(:),thkerr(:),tzerr(:),szerr(:) real(r4),allocatable :: accept_ptemp(:,:),accept_saln(:,:) real(r4),allocatable :: accept_ctemp(:,:),accept_csaln(:,:),accept_ctstd(:,:),accept_csstd(:,:) real(r4),allocatable :: accept_lon(:),accept_lat(:),accept_ptime(:) character(:),allocatable :: var character(cl) :: argoFileNames(3) integer(i4) :: nargoFiles integer(i4) :: total_files_processed ! reject counts integer(i4) :: ndate_time_qc_fail,ninv_qc_fail,nclim_qc_fail real(r4) :: ipr,jpr date_offset=ndays(01,01,1950,12,31,1900) sigver=4 !default total_files_processed=0 ! analysis date !call getarg(1,cdtg) call CDATE2WNDAY(WDAY,cdtg) call read_profile_params() call initialize_analysis_grid() call read_obs_location() call read_obs_window() call initialize_gdem() gdem_data_location=trim(clim_data_location) ! set target z levels tzlevels=[0.,4., 8.,10.,15.,20.,25.,30.,35.,40.,45.,50.,60.,70.,80.,90.,100.,& 110.,120.,130.,140.,150.,160.,170.,180.,190., 200.,210., 220.,230.,& 240.,250.,260.,270., 280.,290.,300.,325.,350.,375.,400.,500.,600.,700.,& 800.,900.,1000.,1200.,1400.,1600.,1800.,2000.] allocate(ftemp1d(ntzlevs),fsaln1d(ntzlevs),ptemp(ntzlevs)) allocate(nlon(nprofmx),nlat(nprofmx),nptemp(nprofmx,ntzlevs),nsaln(nprofmx,ntzlevs),ptime(nprofmx)) nuniq=0 nwithin=0 fcount=0 atlFile=trim(raw_argo_profiles_location)//"atl/"//cdtg(1:8)//"_prof.nc" pacFile=trim(raw_argo_profiles_location)//"pac/"//cdtg(1:8)//"_prof.nc" indFile=trim(raw_argo_profiles_location)//"ind/"//cdtg(1:8)//"_prof.nc" rFile=trim(raw_argo_profiles_location)//"nrt/"//"R"//cdtg(1:8)//"_prof.nc" nargoFiles=0 ! first try delayed time files INQUIRE(FILE=atlFile, EXIST=laexist) INQUIRE(FILE=pacFile, EXIST=lpexist) INQUIRE(FILE=indFile, EXIST=liexist) INQUIRE(FILE=rFile, EXIST=lrexist) if(laexist) then nargoFiles=nargoFiles+1 argoFileNames(nargoFiles)=atlFile endif if(lpexist) then nargoFiles=nargoFiles+1 argoFileNames(nargoFiles)=pacFile endif if(liexist) then nargoFiles=nargoFiles+1 argoFileNames(nargoFiles)=indFile endif dt_argo = laexist .or. lpexist .or. liexist rt_argo = lrexist .and. (.not. dt_argo) if (dt_argo) then print *, "delayed time files available" print *, " will use the foll files: " do i=1,nargoFiles print *, trim(argoFileNames(i)) enddo else if (rt_argo) then print *, "delayed time files not available but real-time file present" print *, "will use the foll: " nargofiles=1 argoFileNames(1)=rFile print *, trim(argoFileNames(1)) else print *, "no argo files" stop endif do ii=1,nargoFiles argoFile=argoFilenames(ii) !wsday=wday+real(ii,r4) !call wday2cdate(wsday,cdtg) !argoFile=trim(argo_profiles_location)//"R"//cdtg(1:8)//"_prof.nc" !print *, trim(argoFile) !INQUIRE(FILE=argoFile, EXIST=lexist) !if(lexist) then ! open and read info from the argo file for the day ! QC flags are in strage character arrays ! read in qc flags and only select data flagged with "1" call nciopn(argoFile,ncid) call ncioin(argoFile,ncid,"N_PROF",nprof) call ncioin(argoFile,ncid,"N_LEVELS",nlvl) if(allocated(lon)) deallocate(lon,lat,juldr8,juld_qc,pos_qc) allocate(lon(nprof),lat(nprof),juldr8(nprof)) allocate(character(len=nprof) :: juld_qc) allocate(character(len=nprof) :: pos_qc) call nciorv(argoFile,ncid,"LONGITUDE",lon) call nciorv(argoFile,ncid,"LATITUDE",lat) call nciorv(argoFile,ncid,"JULD",juldr8) call nciorv(argoFile,ncid,"JULD_QC",juld_qc) call nciorv(argoFile,ncid,"POSITION_QC",pos_qc) ! fcount is the final number of profiles retained ! it is updated and used for file names !fcount=0 do i=1,nprof total_files_processed=total_files_processed+1 if(allocated(pres)) deallocate(pres,temp,saln,pres_qc,temp_qc,saln_qc,valid) allocate(pres(nlvl,1),temp(nlvl,1),saln(nlvl,1),valid(nlvl)) allocate(character(len=nlvl) :: pres_qc) allocate(character(len=nlvl) :: temp_qc) allocate(character(len=nlvl) :: saln_qc) valid=0 ncount=0 if(rt_argo) then call nciorv(argoFile,ncid,"PRES",pres,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"PRES_QC",pres_qc,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"TEMP",temp,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"TEMP_QC",temp_qc,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"PSAL",saln,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"PSAL_QC",saln_qc,[1,i],[nlvl,1]) else if(dt_argo) then call nciorv(argoFile,ncid,"PRES_ADJUSTED",pres,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"PRES_ADJUSTED_QC",pres_qc,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"TEMP_ADJUSTED",temp,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"TEMP_ADJUSTED_QC",temp_qc,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"PSAL_ADJUSTED",saln,[1,i],[nlvl,1]) call nciorv(argoFile,ncid,"PSAL_ADJUSTED_QC",saln_qc,[1,i],[nlvl,1]) endif ! check the QC flags and flag invalid profiles do k=1,nlvl if(pres_qc(k:k)=="1" .and. temp_qc(k:k)=="1" .and. saln_qc(k:k)=="1" ) then valid(k)=1 endif enddo nvalid=sum(valid) ! select files with position qc and time qc with flags 1 and within the ! analysis domain ! add an extra constraint thata minimum of 5 points are there in each ! profile if(nvalid>5.0 .and. lon(i)> lnmn .and. lon(i)ltmn .and. lat(i)1) then prof_valid(i)=0 endif enddo print *, "Date Q/C passed: ",sum(prof_valid) !ndate_time_qc_fail=fcount-sum(prof_valid) ! check for inversion and reject based on drho_tol !Write(*,*) "Checking for inversion and reject based on drho_tol" !write(*,'(a)') "**********************************************" !write(*,'(a)') "Profile No. Depth [m] Inversion magnitude" !write(*,'(a)') "**********************************************" do i=1,fcount if(prof_valid(i)==1) then ninv=0 do k=1,ntzlevs prof_dens(i,k)= sig(nptemp(i,k),nsaln(i,k),sigver) enddo do k=2,ntzlevs if(nptemp(i,k)/=-999. .and. (prof_dens(i,k)-prof_dens(i,k-1))<-drho_tol) then ninv=ninv+1 prof_valid(i)=0 !write(*,'(I5,t15,F10.4,t30,f10.4)') i,tzlevels(k),(prof_dens(i,k)-prof_dens(i,k-1)) endif enddo endif enddo ninv_qc_fail=fcount-sum(prof_valid)-ndate_time_qc_fail print *, "Inversion Q/C passed: ",sum(prof_valid) allocate(atclim(fcount,ntzlevs)) allocate(asclim(fcount,ntzlevs)) allocate(atstd(fcount,ntzlevs)) allocate(asstd(fcount,ntzlevs)) call get_clim_profiles(fcount,nlon,nlat,cdtg,ntzlevs,tzlevels,atclim,asclim,atstd,asstd) !write(*,*) !Write(*,*) "Checking profiles against GDEM-V3 and rejecting values outside 3sig of clim std" !write(*,'(a)') "****************************************************************************************" !write(*,'(a)') "Profile No. Depth [m] [temp-clim] nsig*std [saln-clim] nsig*std" !write(*,'(a)') "****************************************************************************************" !! climatology check prof_kmx=40 do i=1,fcount if(prof_valid(i)==1) then do k=1,ntzlevs if(nsaln(i,k)==-999.0) then kmx=k-1 exit endif enddo prof_kmx(i)=kmx n_outliers=0 do k=1,kmx if((abs(nptemp(i,k)-atclim(i,k))>nsig*atstd(i,k)) .or. (abs(nsaln(i,k)-asclim(i,k))>nsig*asstd(i,k))) then !write(*,'(I4,t15,f12.4,t30,f12.4,t45,f12.4,t60,f12.4,t75,f12.4)') i,tzlevels(k),nptemp(i,k)-atclim(i,k),nsig*atstd(i,k),nsaln(i,k)-asclim(i,k),nsig*asstd(i,k) n_outliers=n_outliers+1 if(n_outliers>frac_outliers*ntzlevs) then prof_valid(i)=0 endif endif enddo endif enddo nclim_qc_fail=fcount-sum(prof_valid)-ndate_time_qc_fail-ninv_qc_fail write(*,*) "Climatological QC Failed: ", nclim_qc_fail print *, "Final no of accepted profiles after climatology check: ",sum(prof_valid) !ndups=nwithin-nuniq !write(*,*) "Date-Time QC Failed: ", ndate_time_qc_fail !write(*,*) "Inversion QC Failed: ", ninv_qc_fail !write(*,*) "Climatological QC Failed: ", nclim_qc_fail !write(*,*) "No of Profiles Processed: ", nwithin !write(*,*) "No of Duplicates Rejected: ", ndups !write(*,*) "Final no of accepted profiles: ",fcount-ndate_time_qc_fail-ninv_qc_fail-nclim_qc_fail !print *, sum(prof_valid) !!transfer the accepted files for output ! calculate the no of accepted profiles ncount=sum(prof_valid) allocate(accept_lon(ncount)) allocate(accept_lat(ncount)) allocate(accept_ptime(ncount)) allocate(accept_ptemp(ncount,ntzlevs)) allocate(accept_saln(ncount,ntzlevs)) allocate(accept_ctemp(ncount,ntzlevs)) allocate(accept_csaln(ncount,ntzlevs)) allocate(accept_ctstd(ncount,ntzlevs)) allocate(accept_csstd(ncount,ntzlevs)) allocate(prof_num(ncount)) vcount=0 do i=1,fcount if(prof_valid(i)==1) then vcount=vcount+1 prof_num(vcount)=vcount accept_lon(vcount)=nlon(i) accept_lat(vcount)=nlat(i) accept_ptime(vcount)=ptime(i) accept_ptemp(vcount,:)=nptemp(i,:) accept_saln(vcount,:)=nsaln(i,:) accept_ctemp(vcount,:)=atclim(i,:) accept_csaln(vcount,:)=asclim(i,:) accept_ctstd(vcount,:)=atstd(i,:) accept_csstd(vcount,:)=asstd(i,:) endif enddo if (vcount >= 1 ) then !! Alex add condition fileName=trim(argo_tsidb_location)//"tsidb_argo_"//cdtg//".nc" call write_argo_tsidb(filename,ncount,ntzlevs,prof_num,tzlevels,accept_lon,accept_lat,accept_ptime,& accept_ptemp,accept_saln,accept_ctemp,accept_csaln,accept_ctstd,accept_csstd) endif end subroutine subroutine argo2tsis_profiles(cdtg) implicit none character(:),allocatable :: argoFile,dbloc,argoFilelist character(10) :: cdtg,cdtgl integer :: ncid,nprof,nlvl,vid,fid real(r4),allocatable :: lon(:),lat(:) integer(i4),parameter :: ntzlevs=52 real(r4) :: tzlevels(ntzlevs) integer :: i,j,ncount,k,kk,stat,ngdem,fcount,ii,jj,kmx,l,nlp character(:),allocatable :: FileName character(:), allocatable :: vname integer(i4) :: vtype integer(i4) :: vdims3(3),vdims2(2),vdim(1) real(r4) :: date_offset,wday,wsday logical(bl) :: lexist integer,parameter :: nprofmx=10000 real(r4),allocatable :: plon(:),plat(:),ptemp(:,:),psaln(:,:),ptime(:),pdens(:,:) integer(i4),allocatable :: indx(:,:),akmx(:,:),prof_kmx(:),prof_valid(:),prof_num(:) real(r4),allocatable :: pl(:),rl(:),sl(:),tl(:),ul(:),vl(:),dp(:),terr(:),serr(:),derr(:),thkerr(:),tzerr(:),szerr(:) real(r4),allocatable :: plt(:),rlt(:),slt(:),tlt(:),iplt(:),tzt1(:),tzt2(:) character(:),allocatable :: var real(r4),parameter :: spval=2.00**100 real(r4) :: rho_min,tsig,q ! error dp or d real(r4),parameter :: sigdmin = 0.005 real(r4),parameter :: sigdpmin = 0.5 real(r4),parameter :: factdp = 0.02 real(r4),parameter :: fact_tserr=0.5 real(r4) :: fact_tserr_tau real(r4),parameter :: argo_temp_err=0.01, argo_sal_err=0.005 real(r4) :: il,jl !model info character(len=20) :: model,fcst_file,anl_file logical(bl) :: lglb,fcst_file_out character(len=5) :: pfx NAMELIST/model_info/model,sigver,lglb,fcst_file,anl_file,fcst_file_out,pfx ! reject counts integer(i4) :: ndate_time_qc_fail,ninv_qc_fail,nclim_qc_fail real(r4) :: ipr,jpr,min_rho date_offset=ndays(01,01,1950,12,31,1900) ! set target z levels tzlevels=[0.,4., 8.,10.,15.,20.,25.,30.,35.,40.,45.,50.,60.,70.,80.,90.,100.,& 110.,120.,130.,140.,150.,160.,170.,180.,190., 200.,210., 220.,230.,& 240.,250.,260.,270., 280.,290.,300.,325.,350.,375.,400.,500.,600.,700.,& 800.,900.,1000.,1200.,1400.,1600.,1800.,2000.] !call getarg(1,cdtg) call CDATE2WNDAY(WDAY,cdtg) call read_profile_params() call initialize_analysis_grid() call read_obs_location() call read_obs_window() call initialize_pobs_vector() OPEN(21,FILE='tsis.nlist') READ(21,NML=model_info) CLOSE(21) ncount=0 !1st pass to get the no of profiles do ii=twin_start,twin_end wsday=wday+real(ii,r4) call wday2cdate(wsday,cdtgl) cdtgl(9:10)="00" argoFile=trim(argo_tsidb_location)//"tsidb_argo_"//cdtgl//".nc" INQUIRE(FILE=argoFile, EXIST=lexist) if(lexist) then call nciopn(argoFile,fid) call ncioin(argoFile,fid,"nprof",nprof) ncount=ncount+nprof call nciocl(argoFile,fid) endif print *, trim(argoFile),nprof,ncount enddo if(ncount>1) then allocate(indx(nx,ny)) allocate(ptime(ncount)) allocate(plon(ncount)) allocate(plat(ncount)) allocate(ptemp(ncount,ntzlevs)) allocate(psaln(ncount,ntzlevs)) allocate(pdens(ncount,ntzlevs)) allocate(tzerr(ntzlevs)) allocate(szerr(ntzlevs)) allocate(tzt1(ntzlevs)) allocate(tzt2(ntzlevs)) allocate(pl(nz+1),rl(nz),sl(nz),tl(nz),ul(nz),vl(nz),dp(nz),derr(nz),thkerr(nz),terr(nz),serr(nz)) allocate(plt(nz+1),rlt(nz),slt(nz),tlt(nz),iplt(nz)) ! read the data in the second pass ncount=0 do ii=twin_start,twin_end wsday=wday+real(ii,r4) call wday2cdate(wsday,cdtgl) cdtgl(9:10)="00" argoFile=trim(argo_tsidb_location)//"tsidb_argo_"//cdtgl//".nc" INQUIRE(FILE=argoFile, EXIST=lexist) if(lexist) then call nciopn(argoFile,fid) call ncioin(argoFile,fid,"nprof",nprof) call nciorv(argoFile,fid,"lon",plon(ncount+1:ncount+nprof),[1],[nprof]) call nciorv(argoFile,fid,"lat",plat(ncount+1:ncount+nprof),[1],[nprof]) call nciorv(argoFile,fid,"z",tzlevels) call nciorv(argoFile,fid,"time",ptime(ncount+1:ncount+nprof),[1],[nprof]) call nciorv(argoFile,fid,"temp",ptemp(ncount+1:ncount+nprof,:),[1,1],[nprof,ntzlevs]) call nciorv(argoFile,fid,"saln",psaln(ncount+1:ncount+nprof,:),[1,1],[nprof,ntzlevs]) ncount=ncount+nprof call nciocl(argoFile,fid) endif ! print *, trim(argoFile),nprof,ncount enddo open(unit=10,file="argo.txt",status="unknown") do i=1,ncount write(10,'(2F12.4)') plon(i),plat(i) enddo close(10) fcount=ncount ncount=0 do i=1,fcount if (plon(i) > lnmn .and. plon(i)ltmn .and.plat(i)