!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module estuary_mod ! !DESCRIPTION: ! This module contains the virtual salt flux version of the ! estuary parameterization ! ! !USES: use POP_KindsMod use POP_IOUnitsMod use kinds_mod use domain_size use domain use constants use io use io_types use grid use forcing_fields,only:ROFF_F use time_management use prognostic,only:TRACER,curtime use global_reductions, only: global_sum_prod use tavg, only: define_tavg_field, accumulate_tavg_field, accumulate_tavg_now implicit none private save ! !PUBLIC MEMBER FUNCTIONS: public :: init_estuary, & set_estuary_vsf_forcing, & set_estuary_exch_circ, & add_estuary_vsf_tracer_tend private :: estuary_box_model ! !PUBLIC DATA MEMBERS: logical, public :: lestuary_on, lvsf_river, lebm_on real(r8), public :: vsf_river_correction ! !PRIVATE DATA MEMBERS: ! Maximum layer index to consider estuary contributions integer :: kmax_estuary = 15 ! Surface Virtual Salt Flux Associated with Rivers real (r8), dimension(nx_block,ny_block,max_blocks_clinic) ::& FLUX_ROFF_VSF_SRF ! Vertical Salt Flux Across Upper/Lower Layer Interface (From box model) REAL (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: & FLUX_EXCH_INTRF ! Vertical salt fluxes across top and bottom of a cell interface real (r8), dimension(nx_block,ny_block,max_blocks_clinic) ::& FLUX_ROFF_VSF, FLUX_ROFF_VSF_KM1, & FLUX_EXCH, FLUX_EXCH_KM1 ! EBM parameters expanded to 2D fields real (r8), dimension(nx_block,ny_block, max_blocks_clinic) :: & & est_tide_amp,& & est_mouth_width, & & est_mouth_depth, & & est_length_a1,& & est_tidal_pump_a2, & & est_lower_depth_ratio, & & est_thick_upper, & & est_thick_lower ! EBM input and output fields real (r8), dimension(nx_block,ny_block, max_blocks_clinic) :: & & S_lower_layer, S_upper_layer, & & Q_river, Q_upper, Q_lower !----------------------------------------------------------------------- ! ! tavg ids for tavg diagnostics related to EBM ! !----------------------------------------------------------------------- integer (int_kind) :: & tavg_FLUX_ROFF_VSF_SRF, &!Surf virt salt flux from river tavg_FLUX_EXCH_INTRF, &!Vert salt flux across layer interface from EBM tavg_FLUX_ROFF_VSF, &!Vert salt flux across cell level interface tavg_FLUX_EXCH !EOC !*********************************************************************** contains !*********************************************************************** subroutine init_estuary !----------------------------------------------------------------------- ! ! initialize estuary parameterization ! !-----------------------------------------------------------------------0 implicit none integer (int_kind) :: & i,j,k,n, &! level index nu, &! unit for input dataset nml_error, &! namelist error flag bid, &! local block address for this block iblock !---- Varaibles for global EBM parameter input file type(datafile) :: ebm_param_file type(io_field_desc) :: io_tide_amp, io_W_h, io_H, io_a1, io_a2, & & io_h0, io_h_upper, io_h_lower type(io_dim) :: i_dim, j_dim character(char_len) :: io_field_name integer(i4) :: io_field_loc, io_field_type real (r8), allocatable, dimension(:,:,:) :: TEMP_DATA ! temporary data array !Set parameters to the same value for all estuaries (simple test) !Eventually replace this with an input file real(r8) :: & tide_amp = 1.0, & !Tidal amplitude for estuary (m/s), used only if a2>0 W_h = 2000.0, & !Estuary width (m) H = 10.0, & !Estuary mouth depth (m) a1 = 0.876, & !Coefficient for estuary length (optimized for Geyer 2010) a2 = 0.0, & !Coefficient for tidal pumping (a2=0 for no tidal pumping) h0 = 0.5 , & !Ratio of lower-layer depth to H at estuary mouth h_upper = 10.0, & ! Thickness of upper layer of exchange flow h_lower = 10.0 ! Thickness of lower layer of exchange flow character (char_len) :: & estuary_option, & ! estuary option (on or off) estuary_type, & ! estuary type (vsf or vsf_ebm) ebm_param_option, & ! estuary box model parameter type (internal or file) ebm_param_filename, & ! estuary box model parameter file ebm_param_file_fmt ! estuary box model parameter file format (nc,bin) namelist /estuary_nml/ & estuary_option, & ! 'on', 'off' estuary_type, & ! 'vsf', 'vsf_ebm' ebm_param_option, & ! 'internal', 'file' ebm_param_filename, & ! full pathname for input file ebm_param_file_fmt, & ! 'nc,bin' tide_amp, W_h, H, a1, a2, h0, h_upper, h_lower ! EBM parameters !----------------------------------------------------------------------- ! Everything off by default estuary_option = 'off' estuary_type = 'vsf' ebm_param_option = 'internal' ebm_param_filename = 'unknown-file' ebm_param_file_fmt = 'nc' ! Look for namelist block if (my_task == master_task) then open (nml_in, file=nml_filename, status='old',iostat=nml_error) if (nml_error /= 0) then nml_error = -1 else nml_error = 1 endif do while (nml_error > 0) read(nml_in, nml=estuary_nml,iostat=nml_error) end do if (nml_error == 0) close(nml_in) endif call broadcast_scalar(nml_error, master_task) if (nml_error /= 0) then call exit_POP(sigAbort,'ERROR reading estuary_nml') endif ! Echo input parameters if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,ndelim_fmt) write(stdout,blank_fmt) write(stdout,*) ' Estuary VSF parameterization:' write(stdout,blank_fmt) write(stdout,*) ' Estuary VSF parametierzation namelist ' write(stdout,blank_fmt) write(stdout, estuary_nml) write(stdout,blank_fmt) endif call broadcast_scalar(estuary_option, master_task) call broadcast_scalar(estuary_type, master_task) call broadcast_scalar(ebm_param_option, master_task) call broadcast_scalar(ebm_param_filename,master_task) call broadcast_scalar(ebm_param_file_fmt,master_task) call broadcast_scalar(tide_amp, master_task) call broadcast_scalar(W_h, master_task) call broadcast_scalar(H, master_task) call broadcast_scalar(a1, master_task) call broadcast_scalar(a2, master_task) call broadcast_scalar(h0, master_task) call broadcast_scalar(h_upper, master_task) call broadcast_scalar(h_lower, master_task) if (estuary_option .EQ. 'on' ) then lestuary_on = .true. elseif (estuary_option .EQ. 'off' ) then lestuary_on = .false. else call exit_POP(sigAbort,'ERROR estuary_option unknown') endif if ( .NOT. lestuary_on ) return if (estuary_type == 'vsf' ) then lvsf_river = .true. lebm_on = .false. elseif ( estuary_type == 'vsf_ebm' ) then lvsf_river = .true. lebm_on = .true. else call exit_POP(sigAbort,'ERROR estuary_type unknown') endif if (ebm_param_option == 'file') then allocate( TEMP_DATA(nx_block,ny_block,max_blocks_clinic)) ! Set up the input file ebm_param_file = construct_file(ebm_param_file_fmt,& & full_name=trim(ebm_param_filename),& & record_length=rec_type_dbl,& & recl_words=nx_global*ny_global) call data_set(ebm_param_file,'open_read') i_dim = construct_io_dim('nx',nx_global) j_dim = construct_io_dim('ny',ny_global) io_field_loc = field_loc_center io_field_type = field_type_scalar ! Grab the data fields and copy them into the internal 2D arrays (est_* ) TEMP_DATA = c0 io_tide_amp = construct_io_field('tide_amp',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_tide_amp) call data_set(ebm_param_file,'read',io_tide_amp) call destroy_io_field(io_tide_amp) est_tide_amp = TEMP_DATA TEMP_DATA = c0 io_W_h = construct_io_field('W_h',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_W_h) call data_set(ebm_param_file,'read',io_W_h) call destroy_io_field(io_W_h) est_mouth_width = TEMP_DATA TEMP_DATA = c0 io_H = construct_io_field('H',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_H) call data_set(ebm_param_file,'read',io_H) call destroy_io_field(io_H) est_mouth_depth = TEMP_DATA TEMP_DATA = c0 io_a1 = construct_io_field('a1',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_a1) call data_set(ebm_param_file,'read',io_a1) call destroy_io_field(io_a1) est_length_a1 = TEMP_DATA TEMP_DATA = c0 io_a2 = construct_io_field('a2',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_a2) call data_set(ebm_param_file,'read',io_a2) call destroy_io_field(io_a2) est_tidal_pump_a2 = TEMP_DATA TEMP_DATA = c0 io_h0 = construct_io_field('h0',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_h0) call data_set(ebm_param_file,'read',io_h0) call destroy_io_field(io_h0) est_lower_depth_ratio = TEMP_DATA TEMP_DATA = c0 io_h_upper = construct_io_field('h_upper',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_h_upper) call data_set(ebm_param_file,'read',io_h_upper) call destroy_io_field(io_h_upper) est_thick_upper = TEMP_DATA TEMP_DATA = c0 io_h_lower = construct_io_field('h_lower',dim1=i_dim,dim2=j_dim,& & field_loc=io_field_loc,field_type=io_field_type,& & d2d_array=TEMP_DATA(:,:,:)) call data_set(ebm_param_file,'define',io_h_lower) call data_set(ebm_param_file,'read',io_h_lower) call destroy_io_field(io_h_lower) est_thick_lower = TEMP_DATA !!! call exit_POP(sigAbort,'ERROR not ready for estuary file input') elseif (ebm_param_option == 'internal' ) then est_tide_amp = tide_amp est_mouth_width = W_h est_mouth_depth = H est_length_a1 = a1 est_tidal_pump_a2 = a2 est_lower_depth_ratio = h0 est_thick_upper = h_upper*100. ! m -> cm est_thick_lower = h_lower*100. ! m -> cm else call exit_POP(sigAbort,'ERROR ebm_param_option unknown') endif !----------------------------------------------------------------------- ! ! define tavg fields related to EBM ! !----------------------------------------------------------------------- call define_tavg_field(tavg_FLUX_ROFF_VSF_SRF,'FLUX_ROFF_VSF_SRF',2, & long_name='Surface Virtual Salt Flux Associated with Rivers',& units='g/kg*cm/s', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_FLUX_EXCH_INTRF,'FLUX_EXCH_INTRF',2, & long_name='Vertical Salt Flux Across Upper/Lower Layer Interface (From EBM)',& units='g/kg*cm/s', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_FLUX_ROFF_VSF,'FLUX_ROFF_VSF',3, & long_name='Vertical salt fluxes across the cell interface from ROFF', & units='g/kg*cm/s', grid_loc='3112', & coordinates='TLONG TLAT z_w time') call define_tavg_field(tavg_FLUX_EXCH,'FLUX_EXCH',3, & long_name='Vertical salt fluxes across the cell interface from exchange circ', & units='g/kg*cm/s', grid_loc='3112', & coordinates='TLONG TLAT z_w time') !----------------------------------------------------------------------- return end subroutine init_estuary !*********************************************************************** !BOP ! !IROUTINE: set_estuary_vsf_forcing ! !INTERFACE: subroutine set_estuary_vsf_forcing ! !DESCRIPTION: ! This routine calucates the salinity flux through the sea surface ! for the virtual salt flux forcing option. It uses the local SSS. ! ! Need to add code to compute correction for global conservation ! ! !REVISION HISTORY: ! same as module ! !INPUT/OUTPUT PARAMETERS: !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: iblock real(r8) :: gsum_vsf_sref, gsum_vsf_sloc !----------------------------------------------------------------------- ! Compute surface virtual salt flux using local salinity do iblock=1,nblocks_clinic FLUX_ROFF_VSF_SRF(:,:,iblock) = -fwmass_to_fwflux*ROFF_F(:,:,iblock)* & & TRACER(:,:,1,2,curtime,iblock)* & & RCALCT(:,:,iblock) ! !QS: diag output ! IF (my_task.EQ.47) THEN ! WRITE(stdout,71)iblock,ROFF_F(34,30,iblock), & ! & TRACER(34,30,1,2,curtime,iblock), & ! & FLUX_ROFF_VSF_SRF(34,30,iblock) !71 FORMAT('QS: tot_ROFF_VSF (iblock,ROFF_F,SSS,ROFF_VSF_SRF): ', & ! & 1(i3,1x),3(1pe11.4,1x)) ! END IF enddo ! Compute correction for global conservation gsum_vsf_sref = global_sum_prod(ROFF_F,TAREA,distrb_clinic,& & field_loc_center, RCALCT)*salinity_factor gsum_vsf_sloc = global_sum_prod(FLUX_ROFF_VSF_SRF,TAREA,distrb_clinic,& & field_loc_center, RCALCT) vsf_river_correction = (gsum_vsf_sref - gsum_vsf_sloc)/area_t !!$ if ( my_task == master_task) then !!$ write(stdout,*) ' ### ESTUARY FORCING ###' !!$ write(stdout,*) ' vsf_sref = ', gsum_vsf_sref !!$ write(stdout,*) ' vsf_sloc = ', gsum_vsf_sloc !!$ write(stdout,*) ' vsf correction = ', vsf_river_correction !!$ endif return end subroutine set_estuary_vsf_forcing !*********************************************************************** !BOP ! !IROUTINE: set_estuary_exch_circ ! !INTERFACE: subroutine set_estuary_exch_circ ! !DESCRIPTION: ! This routine calculates the salinity flux through the interface ! separating the lower and upper layers of the exchange circulation ! ! !REVISION HISTORY: ! same as module ! !INPUT/OUTPUT PARAMETERS: ! !INPUT PARAMETERS: !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- real(r8) :: z_layer_top, z_layer_bot ! depths of top and bottom of lower layer real(r8) :: z_cell_top, z_cell_bot ! depths of top and bottom of grid cell real (r8) :: wgt integer :: i,j,k,iblock !----------------------------------------------------------------------- S_lower_layer = c0 Q_river = fwmass_to_fwflux*ROFF_F*TAREA*1.0e-6 ! kg/m^2/s -> m^3/s FLUX_EXCH_INTRF = c0 do iblock=1,nblocks_clinic do j=1,ny_block do i=1,nx_block if ( est_thick_lower(i,j,iblock) == c0 .OR. & & Q_river(i,j,iblock) == 0 ) cycle ! no estuary at this point z_layer_top = est_thick_upper(i,j,iblock) z_layer_bot = est_thick_upper(i,j,iblock)+est_thick_lower(i,j,iblock) ! Compute lower layer averaged salinity do k=1,kmax_estuary if ( k == 1 ) then z_cell_top = 0. else z_cell_top = zw(k-1) endif z_cell_bot = zw(k) if ( z_cell_bot < z_layer_top ) cycle ! still above lower layer if ( z_cell_top > z_layer_bot ) exit ! now bleow lower layer ! Fraction of the lower layer in this cell wgt = (min(z_cell_bot,z_layer_bot) - & & max(z_cell_top,z_layer_top))/est_thick_lower(i,j,iblock) S_lower_layer(i,j,iblock) = S_lower_layer(i,j,iblock) + & wgt*TRACER(i,j,k,2,curtime,iblock)*c1000 enddo ! Run the estuary box model to get the exchange flow volume flux call estuary_box_model(Q_river(i,j,iblock),& & est_tide_amp(i,j,iblock),S_lower_layer(i,j,iblock),& & est_mouth_width(i,j,iblock),est_mouth_depth(i,j,iblock),& & est_length_a1(i,j,iblock),est_tidal_pump_a2(i,j,iblock),& & est_lower_depth_ratio(i,j,iblock),& & Q_upper(i,j,iblock),Q_lower(i,j,iblock),S_upper_layer(i,j,iblock)) ! Translate the exchange flow volume flux into a vertical salt flux ! NB Q_lower is negative FLUX_EXCH_INTRF(i,j,iblock) = -Q_lower(i,j,iblock)*(1.0e6)* & ! m^3/s -> cm^3/s & ( S_lower_layer(i,j,iblock)*(1.0e-3) & ! psu -> msu & -S_upper_layer(i,j,iblock)*(1.0e-3))/ & ! psu -> msu & TAREA(i,j,iblock)*RCALCT(i,j,iblock) ! ! QS: diag. output ! IF (my_task.EQ.47.AND.i.EQ.34.AND.j.EQ.30) THEN ! WRITE(stdout,81)i,j, Q_river(i,j,iblock), & ! & Q_upper(i,j,iblock),S_upper_layer(i,j,iblock), & ! & Q_lower(i,j,iblock),S_lower_layer(i,j,iblock), & ! & Q_river(i,j,iblock)-Q_upper(i,j,iblock)-Q_lower(i,j,iblock), & ! & Q_upper(i,j,iblock)*S_upper_layer(i,j,iblock) & ! & +Q_lower(i,j,iblock)*S_lower_layer(i,j,iblock) !81 FORMAT('QS: POP<--EBM (i,j,Qr,Qu,Su,Ql,Sl,totQ,totSF)'2(i3,1x),7(1pe11.4,1x)) ! WRITE(stdout,82)TAREA(i,j,iblock),FLUX_EXCH_INTRF(i,j,iblock) !82 FORMAT('QS: tot_EXCH_VSF (Tarea,EXCH_VSF_INIRF): ',2(1pe11.4,1x)) ! END IF enddo enddo enddo return end subroutine set_estuary_exch_circ !*********************************************************************** !BOP ! !IROUTINE: add_estuary_vsf_tracer_tend ! !INTERFACE: subroutine add_estuary_vsf_tracer_tend(T_SOURCE, k, this_block) ! !DESCRIPTION: ! This routine calculates the salinity flux through the bottom of a ! tracer cell due to the vertically distributed surface virtual ! salt flux and due to the estuary exchange flow. The salinity tendency ! is computed as the divergence of the fluxes across the cell. ! Only salinity is affected. Transport of other tracers by the exchange ! flow could be easily added later. ! ! This routine must be called in level-order k=1,2,3... ! ! !REVISION HISTORY: ! same as module ! !INPUT/OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: & T_SOURCE ! source terms for all tracers (to avoid copies) ! contribution to tendency only added to salinity ! (tracer index 2) for now. ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & k ! vertical level index type (block), intent(in) :: & this_block ! block info for this block !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & i,j,bid real (r8), dimension(nx_block,ny_block) :: & WORK1 ! temporary work space !----------------------------------------------------------------------- !EOC ! Skip all of this if we are beyond the maximum allowable depth ! for estuary fluxes if ( k > kmax_estuary ) then return endif bid = this_block%local_id ! Set fluxes at top of top cell if ( k == 1 ) then FLUX_ROFF_VSF_KM1(:,:,bid) = FLUX_ROFF_VSF_SRF(:,:,bid) ! Virtual salt flux from runoff at sea surface FLUX_EXCH_KM1(:,:,bid) = 0 ! Vertical salt flux across interface between endif ! Set virtual salt flux at bottom of cell ! NB zw(k) always > 0, so should not get divide by zero if est_thick_upper = 0 do j=1,ny_block do i=1,nx_block if ( zw(k) <= est_thick_upper(i,j,bid) ) then FLUX_ROFF_VSF(i,j,bid) = FLUX_ROFF_VSF_SRF(i,j,bid)*& & (c1 - zw(k)/est_thick_upper(i,j,bid)) else FLUX_ROFF_VSF(i,j,bid) = c0 endif enddo enddo ! Set exchange flow flux at bottom of cell do j=1,ny_block do i=1,nx_block if ( est_thick_lower(i,j,bid) == c0 ) cycle ! Guard against divide by zero if ( zw(k) <= est_thick_upper(i,j,bid) ) then FLUX_EXCH(i,j,bid) = FLUX_EXCH_INTRF(i,j,bid)*zw(k)/est_thick_upper(i,j,bid) elseif ( zw(k) <= est_thick_upper(i,j,bid)+est_thick_lower(i,j,bid) ) then FLUX_EXCH(i,j,bid) = FLUX_EXCH_INTRF(i,j,bid)*& & (est_thick_upper(i,j,bid)+est_thick_lower(i,j,bid)-zw(k))/est_thick_lower(i,j,bid) else FLUX_EXCH(i,j,bid) = c0 endif enddo enddo ! Make sure there is no flux through the sea floor where(k >= KMT(:,:,bid)) FLUX_ROFF_VSF(:,:,bid) = c0 FLUX_EXCH(:,:,bid) = c0 endwhere ! Compute contribution to tendency from divergence of flux ! Surface Virtual salt flux contribution if ( partial_bottom_cells ) then WORK1 = RCALCT(:,:,bid)*(FLUX_ROFF_VSF_KM1(:,:,bid) - FLUX_ROFF_VSF(:,:,bid))/DZT(:,:,k,bid) else WORK1 = RCALCT(:,:,bid)*(FLUX_ROFF_VSF_KM1(:,:,bid) - FLUX_ROFF_VSF(:,:,bid))*dzr(k) endif T_SOURCE(:,:,2) = T_SOURCE(:,:,2) + WORK1 ! ! QS diag output to check ROFF VSF divergence ! IF (my_task.EQ.47.AND.k.LE.3) THEN ! WRITE (stdout,95)k,FLUX_ROFF_VSF_KM1(34,30,bid),FLUX_ROFF_VSF(34,30,bid),c1/dzr(k) !95 FORMAT('QS: ROFF_VSF at (k,surface_T,bottom_T,dz):',1(i3,1x),3(1pe11.4,1x)) ! WRITE (stdout,91)k,T_SOURCE(34,30,2),WORK1(34,30) !91 FORMAT('QS: ROFF_VSF at Amazon (k,Tsource,WORK1): ',1(i3,1x),2(1pe11.4,1x)) ! ! ELSE IF (my_task.EQ.78) THEN ! WRITE (stdout,92)k,FLUX_ROFF_VSF_KM1(6,21,bid),FLUX_ROFF_VSF(6,21,bid),WORK1(6,21) !92 FORMAT('QS: ROFF_VSF at Columbia (surf,bott,dF/dz): ',1(i3,1x),3(1pe11.4,1x)) ! END IF ! Exchange flow contribution ! QS: CODE CHANGE!!!! flip the sign of WORK1 if ( partial_bottom_cells ) then WORK1 = RCALCT(:,:,bid)*(FLUX_EXCH(:,:,bid) - FLUX_EXCH_KM1(:,:,bid))/DZT(:,:,k,bid) else WORK1 = RCALCT(:,:,bid)*(FLUX_EXCH(:,:,bid) - FLUX_EXCH_KM1(:,:,bid))*dzr(k) endif T_SOURCE(:,:,2) = T_SOURCE(:,:,2) + WORK1 ! ! QS diag output to check ROFF VSF divergence ! IF (my_task.EQ.47.AND.k.LE.3) THEN ! WRITE (stdout,96)k,FLUX_EXCH_KM1(34,30,bid),FLUX_EXCH(34,30,bid),c1/dzr(k) !96 FORMAT('QS: EXCH_VSF at (k,surface_T,bottom_T,dz):',1(i3,1x),3(1pe11.4,1x)) ! WRITE (stdout,93)k,T_SOURCE(34,30,2),WORK1(34,30) !93 FORMAT('QS: EXCH_VSF at Amazon (k,Tsource,WORK1): ',1(i3,1x),2(1pe11.4,1x)) ! ELSE IF (my_task.EQ.78) THEN ! WRITE (stdout,94)k,FLUX_EXCH_KM1(6,21,bid),FLUX_EXCH(6,21,bid),WORK1(6,21) !94 FORMAT('QS: EXCH_VSF at Columbia (surf,bott,dF/dz): ',1(i3,1x),3(1pe11.4,1x)) ! END IF !----------------------------------------------------------------------- ! ! compute diagnostics if necessary. ! !----------------------------------------------------------------------- if (k == 1) then call accumulate_tavg_field(FLUX_ROFF_VSF_SRF(:,:,bid), & tavg_FLUX_ROFF_VSF_SRF,bid,1) call accumulate_tavg_field(FLUX_EXCH_INTRF(:,:,bid), & tavg_FLUX_EXCH_INTRF,bid,1) endif ! QS: keep the vsf value at the top of the cell call accumulate_tavg_field(FLUX_ROFF_VSF_KM1(:,:,bid),tavg_FLUX_ROFF_VSF,bid,k) call accumulate_tavg_field(FLUX_EXCH_KM1(:,:,bid),tavg_FLUX_EXCH,bid,k) !----------------------------------------------------------------------- ! QS: then move the bottom vsf to the top of the cell ! Move lower face flux to upper face for next call FLUX_ROFF_VSF_KM1(:,:,bid) = FLUX_ROFF_VSF(:,:,bid) FLUX_EXCH_KM1(:,:,bid) = FLUX_EXCH(:,:,bid) return end subroutine add_estuary_vsf_tracer_tend !*********************************************************************** !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !No modifications have been made for this application (MMWhitney 7/22/15) subroutine estuary_box_model (Q_r,tide_amp,S_l,W_h,H,a1,a2,h0, & Q_u,Q_l,S_u) ! !DESCRIPTION: ! This subroutine is a test version of estaury Box Model v2_4_FORTRA ! EBM v2.4 FORTRAN includes river discharge, tides. ! EBM v2.4 FORTRAN does not include wind mixing, thermal expension ! and precipitation/evaporation. !----------------------------------------------------------------------- ! Estuary Box: ! ! ^ z ! |____________________ _______ ! | | | ! Q_u,S_u <---+-- upper layer <-+-- Q_r | ! |--------------------| -+- | H ! Q_l,S_l -+-> lower layer | | h_l | ! ___|____________________| _|_____| ! x 0 -LE ! ! NOTE: the negative lower layer volume flux leaves POP, positive one ! flows into POP. !----------------------------------------------------------------------- ! Governing Equations: ! The Estuary Box Model (EBM) is assumed steady state, flat bottom and ! flat surface with straight channel and rectangular cross-section. ! The EBM is built basic on three global conservations, water mass ! conservation, water volume conservation and potential energy ! conservation. Combine water mass and volume conservation with linear ! equation of state, the salinity flux conservation can be obtained. ! >>> Water mass conservation: ! rho_r*Q_r = rho_u*Q_u + rho_l*Q_l + (rho_l - rho_u)*a2*Q_ut/2 ! >>> Water volume conservation: ! Q_r = Q_u +Q_l ! >>> Potential energy conservation (PEF: potential energy flux): ! PEF_AD + PEF_HD + PEF_VD + PEF_LF = 0 ! --> PEF_AD=1/2*g*rho_l*Q_l*h+1/2*g*(rho_u*Q_u-rho_r*Q_r)*(H+h) ! --> PEF_HD=1/2*g*a2*(rho_l-rho_u)Q_ut*(H+h) ! --> PEF_VD=-1/2*g*KV*a1*LE*(rho_l-rho_u)*(rho_l+rho_u-2*rho_r)/(rho_l-rho_r)*W ! --> PEF_LF=-1/4*g*((rho_u^2+2*rho_l*rho_r-2*rho_u*rho_r)*(H-h)-H*rho_r^2+h*rho_l^2)/(rho_l-rho_r)*Q_l*h ! !----------------------------------------------------------------------- ! QSun, MMWhitney, 2015.04.20 implicit none ! input variables real (r8), intent(in) :: & Q_r, & !River discharge (m**3/s) tide_amp, & !Averaged tidal amplitude at estuary mouth (m) S_l, & !Salinity at estuary lower layer (ppt) W_h, & !Estuary head width (m) H, & !Estuary averaged depth (m) a2, & !A constant (-), of tidal diffusion a1, & !A constant (-), of estuarine mixing length h0 !A constant (-), of ratio of geometry: h_l/H ! output variables real (r8), intent(out) :: & Q_u, & !Upper layer volume flux (m**3/s, positive) Q_l, & !Lower layer volume flux (m**3/s, negative) S_u !Salinity at estuary upper layer (ppt) real (r8) :: & ERR_EBMv24 !Closing error of EBMv2.4 ! local variables real (r8) :: & rho_r, & !River water density (kg/m**3) rho_l, & !Lower layer inflow density / ocean water density (kg/m**3) rho_u, & !Upper layer inflow density / ocean water density (kg/m**3) u_t, & !Tidal current amplitude near bottom (m/s) u_r, & !Riverine velocity at head of EBM (m/s) u_l, & !Estuarine inflow velocity at mouth of EBM (m/s) u_u, & !Estuarine outflow velocity at mouth of EBM (m/s) u_bar, & !Layer averaged net velocity in EBM (m/s) c_wave, & !Densimetric wave phase speed (m/s) ur0, & !Densimetric riverine Froude number (-) ut0, & !Densimetric tidal current Froude number (-) ul0, & !Densimetric lower layer inflow Froude number (-) uu0, & !Densimetric upper layer outflow Froude number (-) R0, & !Layer densimetric riverine Froude number (-) T0, & !Layer densimetric tidal current Froude number (-) h_l, & !Lower layer water depth of EBM (m) a,b,c,d, & !parameters to solve the cubic equation AD,HD,VD,LF !Different mechanism of salt mass transport in EBM (kg/s) integer (int_kind) :: i,n !Indexes integer (int_kind), dimension(3) :: mask !Arrow real (r8), dimension(3,2) :: & roots !roots of 3rd order polynomial: !Dimension one: real part of the root, !Dimension two: imaginary part of the root. real (r8), parameter :: & g = 9.8_r8, & !Gravity (m/s**2) pi = 3.1415927_r8, & !circumference rho_0 = 1000.0_r8, & !water reference density (km/m**3) beta_S = 7.7e-4_r8, & !Saline contraction coefficient (1/ppt) Sc = 2.2_r8 !Schmidt number (-) ! Executable Statements !----Water density---- rho_r = rho_0 rho_l = rho_0*(1.0_r8+beta_S*S_l) !----River and tide---- u_t = -tide_amp*sqrt(g/H) ! Tidal velocity (toward river) u_r = Q_r/(W_h*H*(1.0_r8-h0)) ! Riverine velocity at head of upper layer c_wave = sqrt(beta_S*S_l*g*H) ! Densimetric wave phase speed (m/s) !----governing equations---- !solve upper outflow: Q_l ur0=u_r/c_wave ut0=u_t/c_wave R0=ur0*(1.0_r8-h0) T0=ut0*(1.0_r8-h0)/pi a =-h0**3.0_r8 b = 2.0_r8*h0**2.0_r8*((2.0_r8-h0)*R0-a2*T0) c = 0.096_r8*a1*h0*(Sc**2.0_r8*R0)**(-1.0_r8/3.0_r8)*R0 & -h0*((2.0_r8-h0)*R0*(R0-2.0_r8*a2*T0)+a2**2.0_r8*T0**2.0_r8) d =-0.048_r8*a1*(Sc**2.0_r8*R0)**(-1.0_r8/3.0_r8) & *R0*(R0-2.0_r8*a2*T0) call cubsolve(b/a,c/a,d/a,roots) !----find the solution of dimensionless lower layer inflow mask=0 n=0 do i=1,3 if (roots(i,1) .LT. 0.0_r8 .AND. roots(i,2) .EQ. 0.0_r8) then mask(i)=1 n=n+1 end if end do if (n .eq. 0) then write(stdout,100)roots(1,1),roots(1,2),roots(2,1),roots(2,2),roots(3,1),roots(3,2) 100 format('EBM error NO solution: ',6(1pe11.4,1x)) ul0=0.0_r8 else if (n .eq. 1) then ul0=SUM(roots(1:3,1)*mask) else if (n .gt. 1) then write(stdout,101)roots(1,1),roots(1,2),roots(2,1),roots(2,2),roots(3,1),roots(3,2) 101 format('EBM error Multiple solution: ',6(1pe11.4,1x)) ul0=0.0_r8 end if !----solving for the upper layer salinity and volume flux at EBM mouth uu0=R0/(1.0_r8-h0)-h0/(1.0_r8-h0)*ul0 S_u=(-S_l*ul0*h0-S_l*a2*T0)/(R0-ul0*h0-a2*T0) Q_l=ul0*h0*H*W_h*c_wave Q_u=uu0*(1.0_r8-h0)*H*W_h*c_wave ! write(stdout,102)Q_r,Q_u,Q_l,S_u,S_l !102 format('EBM--->POP original: ',5(1pe11.4,1x)) !----check if the solution of EBMv24 is colsed u_l=ul0*c_wave u_u=uu0*c_wave u_bar=Q_r/(W_h*H) h_l=H*h0 rho_u=rho_0*(1.0_r8+beta_S*S_u) AD= 0.5_r8*g*rho_l*u_l*h_l**2_r8 & +0.5_r8*g*(rho_u*u_u-rho_r*u_r)*(H**2.0_r8-h_l**2.0_r8) HD=-0.5_r8*a2*g*(rho_l-rho_u)*(H**2.0_r8-h_l**2.0_r8)*u_t/pi VD=-0.5_r8*g*(rho_u-rho_l) & *(rho_l+rho_u-2.0_r8*rho_r)/(rho_l-rho_r) & *0.024_r8*a1*H**2.0_r8 & *(c_wave**4.0_r8/(u_bar*Sc**2.0_r8))**(1.0_r8/3.0_r8) LF=-0.25_r8*g*Q_l/W_h & *( (rho_u**2.0_r8+2.0_r8*rho_l*rho_r-2.0_r8*rho_u*rho_r) & *(H-h_l) & -rho_r**2.0_r8*H+rho_l**2.0_r8*h_l)/(rho_l-rho_r) ERR_EBMv24=AD-HD-VD-LF ! write(stdout,103)AD,HD,VD,LF,ERR_EBMv24 !103 format('Salt budget in EBM: ',5(1pe11.4,1x)) ! effective salinity for POP using, the Q_l is a negative real value S_u=-Q_l*S_l/Q_u ! write(stdout,104)Q_r,Q_u,Q_l,S_u,S_l !104 format('EBM--->POP adjusted: ',5(1pe11.4,1x)) end subroutine estuary_box_model !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !No modifications have been made for this application (MMWhitney 7/22/15) subroutine cubsolve(a,b,c,roots) ! !DESCRIPTION: ! This subroutine is to solve the 3rd order polynomial analytically. ! Source of code: http://www.helioscorner.com ! ! QSun, MMWhitney, 2014.06.18 implicit none ! input variables real (r8), intent(in) :: & a,b,c !parameters to solve the cubic equation ! output variables real (r8), dimension(3,2), intent(out) :: & roots !roots of 3rd order polynomial: !Column one: real part of the root, !Column two: imaginary part of the root. ! local variables real (r8) :: & Q, & !Parameter to solve the cubic equation R, & !Parameter to solve the cubic equation Rsqu, & !Parameter to solve the cubic equation Qcub, & !Parameter to solve the cubic equation SQ, & !Parameter to solve the cubic equation theta, & !Parameter to solve the cubic equation X, & !Parameter to solve the cubic equation Y, & !Parameter to solve the cubic equation XY !Parameter to solve the cubic equation real (r8), parameter :: & pi = 3.1415927_r8 !Circumference ! Executable Statements Q = (a**2.0_r8-3.0_r8*b)/9.0_r8 R = (2.0_r8*a**3.0_r8-9.0_r8*a*b+27.0_r8*c)/54.0_r8 Rsqu = R**2.0_r8 Qcub = Q**3.0_r8 if (Rsqu .LT. Qcub) then !There are three real roots theta=acos(R/sqrt(Qcub)) SQ=sqrt(Q) roots(1,1) = -2.0_r8*SQ*cos(theta/3.0_r8)-a/3.0_r8 roots(2,1) = -2.0_r8*SQ*cos((theta+2.0_r8*pi)/3.0_r8)-a/3.0_r8 roots(3,1) = -2.0_r8*SQ*cos((theta-2.0_r8*pi)/3.0_r8)-a/3.0_r8 roots(1,2) = 0.0_r8 roots(2,2) = 0.0_r8 roots(3,2) = 0.0_r8 return endif !Otherwise, there are one real root and two conjugate complex roots X=-(abs(R)+sqrt(Rsqu-Qcub))**(1.0_r8/3.0_r8) if (R .LT. 0.0_r8) then X = -X endif if (X .EQ. 0.0_r8) then Y = 0.0_r8 else Y = Q/Y endif XY = X+Y roots(1,1) = XY-a/3.0_r8 roots(1,2) = 0.0_r8 roots(2,1) = -0.5_r8*XY-a/3.0_r8 roots(3,1) = -0.5_r8*XY-a/3.0_r8 roots(2,2) = sqrt(3.0_r8)*(X-Y)/2.0_r8 roots(3,2) = -sqrt(3.0_r8)*(X-Y)/2.0_r8 end subroutine cubsolve end module estuary_mod