Skip to content

Commit

Permalink
Added station outputs for all tracers
Browse files Browse the repository at this point in the history
  • Loading branch information
@josephzhang8
josephzhang8 committed May 17, 2024
1 parent c321c67 commit 56b88ae
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 57 deletions.
2 changes: 1 addition & 1 deletion sample_inputs/station.in
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
1 1 1 1 1 1 1 1 1 !on (1)|off(0) flags for elev, air pressure, windx, windy, T, S, u, v, w
1 1 1 1 1 1 1 1 1 !on (1)|off(0) flags for elev,air pressure,windx,windy,T,S,u,v,w,rest of tracers (expanded into subclasses of each module)
4 !# of stations
1 6.5833 54.0000 -1 !Format: station #,x,y,z; if ics=2, x,y are degr in lon/lat. z is z-coord (not distance from surface!). For 3D variables, code will extrapolate above surface/below bottom if necessary
2 7.1583 55.195 -5
Expand Down
34 changes: 2 additions & 32 deletions src/Hydro/schism_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -961,14 +961,12 @@ subroutine schism_init(iorder,indir,iths,ntime)

! Kriging option
! Choice of generalized covariance fucntion
! call get_param('param.in','kr_co',1,kr_co,tmp,stringvalue)
if(kr_co<=0.or.kr_co>4) then
write(errmsg,*)'Wrong kr_co:',kr_co
call parallel_abort(errmsg)
endif

!... Max. for vel. magnitude
! call get_param('param.in','rmaxvel',2,itmp,tmp,stringvalue)
if(rmaxvel<1._rkind) then
write(errmsg,*)'Illegal rmaxvel:',rmaxvel
call parallel_abort(errmsg)
Expand All @@ -978,7 +976,6 @@ subroutine schism_init(iorder,indir,iths,ntime)
rmaxvel2=rmaxvel*1.013_rkind !v-vel

!... min. vel for invoking btrack and for abnormal exit in quicksearch
! call get_param('param.in','velmin_btrack',2,itmp,velmin_btrack,stringvalue)
if(velmin_btrack<=0._rkind) then
write(errmsg,*)'Illegal velmin_btrack:',velmin_btrack
call parallel_abort(errmsg)
Expand All @@ -988,14 +985,12 @@ subroutine schism_init(iorder,indir,iths,ntime)
! to avoid underflow. This should not need to be adjusted
! normally; may need to lower it for some benchmark tests
! Default: btrack_nudge=1.013e-3
! call get_param('param.in','btrack_nudge',2,itmp,btrack_nudge,stringvalue)
if(btrack_nudge<=0._rkind.or.btrack_nudge>0.5_rkind) then
write(errmsg,*)'Illegal btrack_nudge:',btrack_nudge
call parallel_abort(errmsg)
endif

! Test btrack alone (1: rotating Gausshill) - can only be used with pure b-tropic model and pure tri
! call get_param('param.in','ibtrack_test',1,ibtrack_test,tmp,stringvalue)
if(ibtrack_test<0.or.ibtrack_test>1) then
write(errmsg,*)'Illegal ibtrack_test:',ibtrack_test
call parallel_abort(errmsg)
Expand All @@ -1004,7 +999,6 @@ subroutine schism_init(iorder,indir,iths,ntime)
if(ibtrack_test==1.and.lhas_quad) call parallel_abort('INIT: btrack cannot have quads')

! Rouse test
! call get_param('param.in','irouse_test',1,irouse_test,tmp,stringvalue)
if(irouse_test/=0.and.irouse_test/=1) then
write(errmsg,*)'Illegal irouse_test:',irouse_test
call parallel_abort(errmsg)
Expand All @@ -1019,7 +1013,6 @@ subroutine schism_init(iorder,indir,iths,ntime)
endif

!... Inundation algorithm flag (1: better algorithm for fine resolution)
! call get_param('param.in','inunfl',1,inunfl,tmp,stringvalue)
if(inunfl/=0.and.inunfl/=1) then
write(errmsg,*)'Illegal inunfl:',inunfl
call parallel_abort(errmsg)
Expand All @@ -1028,17 +1021,12 @@ subroutine schism_init(iorder,indir,iths,ntime)
!... Numerical shoreline flag (boundary between dry and wet elements)
! 1: the wave forces at the shoreline are set equal to the barotropic gradient, which is calculated on the wet element
! 0: the wave forces at the shoreline are calculated using the hgrad_nodes routine (non physical velocities in very shallow cells)
! call get_param('param.in','shorewafo',1,shorewafo,tmp,stringvalue)
if(shorewafo/=0.and.shorewafo/=1) then
write(errmsg,*)'Illegal shorewafo:',shorewafo
call parallel_abort(errmsg)
endif

! write mode; not used really
!! call get_param('param.in','iwrite',1,iwrite,tmp,stringvalue)

! Elev. i.c. option (elev.ic)
! call get_param('param.in','ic_elev',1,ic_elev,tmp,stringvalue)
if(ic_elev/=0.and.ic_elev/=1) then
write(errmsg,*)'Illegal ic_elev:',ic_elev
call parallel_abort(errmsg)
Expand All @@ -1052,13 +1040,10 @@ subroutine schism_init(iorder,indir,iths,ntime)
endif

! Inverse barometric effects on elev. b.c.
! call get_param('param.in','inv_atm_bnd',1,inv_atm_bnd,tmp,stringvalue)
if(inv_atm_bnd/=0.and.inv_atm_bnd/=1) then
write(errmsg,*)'Illegal inv_atm_bnd:',inv_atm_bnd
call parallel_abort(errmsg)
endif
!Reference atmos. pressure
! call get_param('param.in','prmsl_ref',2,itmp,prmsl_ref,stringvalue)

! Scales for dimensioning in inter-subdomain btrack
! mxnbt=s1_mxnbt*nmm*nvrt is the dimension of btlist (nmm is the max. of all
Expand All @@ -1067,27 +1052,19 @@ subroutine schism_init(iorder,indir,iths,ntime)
! (nbt is the initial # of inter-subdomain trajectories), and
! mnbt*nnbr is the dimension of btrecvq() in routine inter_btrack (nnbr is #
! of nbr processes).
! call get_param('param.in','s1_mxnbt',2,itmp,s1_mxnbt,stringvalue)
! call get_param('param.in','s2_mxnbt',2,itmp,s2_mxnbt,stringvalue)
if(s1_mxnbt<=0._rkind.or.s2_mxnbt<=0._rkind) then
write(errmsg,*)'Illegal s[12]_mxnbt:',s1_mxnbt,s2_mxnbt
call parallel_abort(errmsg)
endif

! Station output option (/=0: need station.in)
! If ics=2, the coord. in station.in must be in lat/lon (degrees)
! call get_param('param.in','iout_sta',1,iout_sta,tmp,stringvalue)

if(iout_sta/=0) then
! call get_param('param.in','nspool_sta',1,nspool_sta,tmp,stringvalue) !output skip
if(nspool_sta<=0) call parallel_abort('Wrong nspool_sta')
if(mod(nhot_write,nspool_sta)/=0) call parallel_abort('mod(nhot_write,nspool_sta)/=0')
!'
endif

!... Read harmonic analysis information (Adapted from ADCIRC)
! call get_param('param.in','iharind',1,iharind,tmp,stringvalue)

!... WWM
! Coupling flag
! 0: decoupled so 2 models will run independently;
Expand All @@ -1099,20 +1076,13 @@ subroutine schism_init(iorder,indir,iths,ntime)
call parallel_abort(errmsg)
endif

!... ramp for the wave forces
! if(nrampwafo/=0.and.nrampwafo/=1) then
! write(errmsg,*)'Unknown nrampwafo',nrampwafo
! call parallel_abort(errmsg)
! endif

! Coupling interval (# of time steps)
if(nstep_wwm<1) then
write(errmsg,*)'Wrong coupling interval:',nstep_wwm
call parallel_abort(errmsg)
endif

! Wave boundary layer option
! call get_param('param.in','iwbl',1,iwbl,tmp,stringvalue)
if(iwbl<0.or.iwbl>2) then
write(errmsg,*)'Wrong iwbl:',iwbl
call parallel_abort(errmsg)
Expand Down Expand Up @@ -3907,15 +3877,15 @@ subroutine schism_init(iorder,indir,iths,ntime)
! Station output option (/=0: need station.in)
! If ics=2, the coord. in station.in must be in lat/lon (degrees)
if(iout_sta/=0) then
nvar_sta=9 !# of output variables
nvar_sta=9+ntracers-2 !# of output variables
if(iorder==0) then
allocate(iof_sta(nvar_sta),stat=istat)
if(istat/=0) call parallel_abort('Sta. allocation failure (1)')
endif

if(myrank==0) then
open(32,file=in_dir(1:len_in_dir)//'station.in',status='old')
! Output variables in order: elev, air pressure, windx, windy, T, S, u, v, w
! Output variables in order: elev, air pressure, windx, windy, T, S, u, v, w, rest of tracers
read(32,*)iof_sta(1:nvar_sta) !on-off flag for each variable
read(32,*)nout_sta
! Following is needed for dimension of nwild2
Expand Down
49 changes: 25 additions & 24 deletions src/Hydro/schism_step.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9974,37 +9974,38 @@ subroutine schism_step(it)

do i=1,nout_sta
ie=iep_sta(i)
if(ie==0) then !not parent
if(ie==0) then !no parent in this rank
iep_flag(i)=0 !for comm. later
sta_out(i,j)=0.d0
sta_out3d(:,i,j)=0.d0
zta_out3d(:,i,j)=0.d0
else !is parent
iep_flag(i)=1
sta_out(i,j)=0.d0 !initialize
select case(j)
case(1) !elev.
swild2(1,1:i34(ie))=eta2(elnode(1:i34(ie),ie))
case(2) !air pressure
swild2(1,1:i34(ie))=pr(elnode(1:i34(ie),ie))
case(3) !wind x
swild2(1,1:i34(ie))=windx(elnode(1:i34(ie),ie))
case(4) !wind y
swild2(1,1:i34(ie))=windy(elnode(1:i34(ie),ie))
case(5) !T
swild2(1:nvrt,1:i34(ie))=tr_nd(1,1:nvrt,elnode(1:i34(ie),ie))
case(6) !S
swild2(1:nvrt,1:i34(ie))=tr_nd(2,1:nvrt,elnode(1:i34(ie),ie))
case(7) !u
if(j==1) then !elev.
swild2(1,1:i34(ie))=eta2(elnode(1:i34(ie),ie))
else if(j==2) then !air pressure
swild2(1,1:i34(ie))=pr(elnode(1:i34(ie),ie))
else if(j==3) then !wind x
swild2(1,1:i34(ie))=windx(elnode(1:i34(ie),ie))
else if(j==4) then !wind y
swild2(1,1:i34(ie))=windy(elnode(1:i34(ie),ie))
else if(j==5) then !T
swild2(1:nvrt,1:i34(ie))=tr_nd(1,1:nvrt,elnode(1:i34(ie),ie))
else if(j==6) then !S
swild2(1:nvrt,1:i34(ie))=tr_nd(2,1:nvrt,elnode(1:i34(ie),ie))
else if(j==7) then !u
!Error: may not be accurate near poles as pframe changes rapidly there
swild2(1:nvrt,1:i34(ie))=uu2(1:nvrt,elnode(1:i34(ie),ie))
case(8) !v
swild2(1:nvrt,1:i34(ie))=vv2(1:nvrt,elnode(1:i34(ie),ie))
case(9) !w
swild2(1:nvrt,1:i34(ie))=ww2(1:nvrt,elnode(1:i34(ie),ie))
case default
call parallel_abort('MAIN: unknown sta. output')
end select
swild2(1:nvrt,1:i34(ie))=uu2(1:nvrt,elnode(1:i34(ie),ie))
else if(j==8) then !v
swild2(1:nvrt,1:i34(ie))=vv2(1:nvrt,elnode(1:i34(ie),ie))
else if(j==9) then !w
swild2(1:nvrt,1:i34(ie))=ww2(1:nvrt,elnode(1:i34(ie),ie))
else if(j<=9+ntracers-2) then
swild2(1:nvrt,1:i34(ie))=tr_nd(j-7,1:nvrt,elnode(1:i34(ie),ie))
else
call parallel_abort('STEP: unknown sta. output')
endif !j

if(j<=4) then !2D var.
sta_out(i,j)=sum(arco_sta(i,1:i34(ie))*swild2(1,1:i34(ie)))
Expand Down Expand Up @@ -10082,7 +10083,7 @@ subroutine schism_step(it)
zta_out3d_gb(:,j,i)=-9999.d0
endif
else
sta_out_gb(j,i)=sta_out_gb(j,i)/nwild2(j)
sta_out_gb(j,i)=sta_out_gb(j,i)/dble(nwild2(j))
if(i>4) then !3D only
sta_out3d_gb(:,j,i)=sta_out3d_gb(:,j,i)/dble(nwild2(j))
zta_out3d_gb(:,j,i)=zta_out3d_gb(:,j,i)/dble(nwild2(j))
Expand Down

0 comments on commit 56b88ae

Please sign in to comment.