subroutine gosbc (is, ie, js, je)
#if defined O_mom || defined O_embm
!=======================================================================
! calculate the average fluxes for next ocean time step
!=======================================================================
implicit none
integer ie, is, je, js, i, j, nc
integer iem1, isp1, jem1, jsp1, k
real f1, f1a, f1l, fh, fs, fwcflx, fwaflx, time
real area, tarea, tsflx, rsocn, tmp, fg, fgs, sulphfac
real area_oa, tarea_oa, oa_add, nn, time_new, T_ph_aoa, T_pc_aoa
real T_period, tarea_oa_t
logical ex
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "calendar.h"
include "csbc.h"
include "coord.h"
include "grdvar.h"
include "tmngr.h"
include "switch.h"
# if defined O_embm
include "cembm.h"
include "atm.h"
# endif
# if defined O_mom
include "mw.h"
# endif
# if defined O_ice
# if defined O_ice_cpts
include "cpts.h"
# endif
include "ice.h"
# endif
# if defined O_mtlm
include "mtlm.h"
# endif
include "levind.h"
# if defined O_fwa
include "fwa.h"
include "cregin.h"
include "tmngr.h"
# endif
# if defined O_sed && !defined O_sed_uncoupled
include "sed.h"
# endif
# if defined O_sulphate_data || defined O_sulphate_data_transient
include "insolation.h"
real tot, twt, cz(1), zero(1)
# endif
# if defined O_embm_read_sflx || defined O_embm_write_sflx
integer ntrec
save ntrec
data ntrec /0/
# endif
isp1 = is + 1
iem1 = ie - 1
jsp1 = js + 1
jem1 = je - 1
# if defined O_mom && defined O_embm
f1 = 1./atatm
fh = 2.389e-8/atatm
fs = -socn/atatm
!-----------------------------------------------------------------------
! calculate average net fluxes. convert heat flux to cal/cm**2/s
! from kW/m**2 and fresh water flux (cm/s) to an apparent salt
! flux (g/cm**2/s) using global ocean average salinity, socn
!-----------------------------------------------------------------------
do j=2,jmtm1
do i=2,imtm1
if (tmsk(i,j) .ge. 0.5) then
# if defined O_plume
# if !defined O_plume_brine
subflux(i,j,1) = flux(i,j,isat)
if (subflux(i,j,1) .gt. 0.) subflux(i,j,1) = 0.
flux(i,j,isat) = flux(i,j,isat) - subflux(i,j,1)
subflux(i,j,1) = fh*subflux(i,j,1)
subflux(i,j,2) = flux(i,j,ishum)
# endif
if (subflux(i,j,2) .gt. 0.) subflux(i,j,2) = 0.
flux(i,j,ishum) = flux(i,j,ishum) - subflux(i,j,2)
subflux(i,j,2) = fs*subflux(i,j,2)
# endif
# if defined O_convect_brine
cba0(i,j) = 0.
do nc=0,ncat
cba(i,j,nc) = f1*cba(i,j,nc)
if (cba(i,j,nc) .gt. 0.) then
cbf(i,j,nc) = fs*cbf(i,j,nc)/cba(i,j,nc)
cba0(i,j) = cba0(i,j) + cba(i,j,nc)
else
cbf(i,j,nc) = 0.
cba(i,j,nc) = 0.
endif
enddo
if (cba0(i,j) .gt. 1.) then
if (cba0(i,j) .gt. 1.000001) then
print*, "==> Warning: ice area > 1: ", cba0(i,j)
endif
cba0(i,j) = 1.
endif
cba0(i,j) = 1. - cba0(i,j)
# endif
sbc(i,j,ihflx) = sbc(i,j,ihflx) + fh*flux(i,j,isat)
! add virtual fluxes of salinity
sbc(i,j,isflx) = sbc(i,j,isflx) + fs*flux(i,j,ishum)
else
# if defined O_plume
subflux(i,j,1) = 0.
subflux(i,j,2) = 0.
# endif
# if defined O_convect_brine
cbf(i,j,:) = 0.
cba(i,j,:) = 0.
cba0(i,j) = 1.
# endif
sbc(i,j,ihflx) = 0.
sbc(i,j,isflx) = 0.
endif
if (umsk(i,j) .ge. 0.5) then
# if defined O_ice_evp || defined O_embm_awind
sbc(i,j,itaux) = f1*flux(i,j,nat+1)
sbc(i,j,itauy) = f1*flux(i,j,nat+2)
# endif
else
sbc(i,j,itaux) = 0.
sbc(i,j,itauy) = 0.
endif
enddo
enddo
call setbcx (sbc(1,1,ihflx), imt, jmt)
call setbcx (sbc(1,1,isflx), imt, jmt)
call setbcx (sbc(1,1,itaux), imt, jmt)
call setbcx (sbc(1,1,itauy), imt, jmt)
# if defined O_convect_brine
do nc=0,ncat
call setbcx (cbf(1,1,nc), imt, jmt)
call setbcx (cba(1,1,nc), imt, jmt)
enddo
call setbcx (cba0, imt, jmt)
# endif
# endif
!-----------------------------------------------------------------------
! update boundary conditions from the land model
! do this now instead of in gasbc so fields can be written out
!-----------------------------------------------------------------------
f1l = 0.
f1a = 0.
# if defined O_embm
if (atatm .ne. 0.) f1a = 1.0/atatm
# endif
# if defined O_mtlm
if (atlnd .ne. 0.) f1l = 1.0/atlnd
do j=2,jmtm1
do i=2,imtm1
if (land_map(i,j) .ne. 0) then
sbc(i,j,iro) = sbc(i,j,iro)*f1l
sbc(i,j,isca) = sbc(i,j,isca)*f1l
sbc(i,j,ievap) = sbc(i,j,ievap)*f1l
sbc(i,j,ilwr) = sbc(i,j,ilwr)*f1l
sbc(i,j,isens) = sbc(i,j,isens)*f1l
# if defined O_carbon
sbc(i,j,inpp) = sbc(i,j,inpp)*f1l
sbc(i,j,isr) = sbc(i,j,isr)*f1l
# endif
else
sbc(i,j,iro) = sbc(i,j,iro)*f1a
sbc(i,j,ievap) = 0.
sbc(i,j,ilwr) = 0.
sbc(i,j,isens) = 0.
# if defined O_mtlm && defined O_carbon
sbc(i,j,inpp) = 0.
sbc(i,j,isr) = 0.
# endif
endif
enddo
enddo
call setbcx (sbc(1,1,isca), imt, jmt)
call setbcx (sbc(1,1,ievap), imt, jmt)
call setbcx (sbc(1,1,ilwr), imt, jmt)
call setbcx (sbc(1,1,isens), imt, jmt)
# if defined O_carbon
call setbcx (sbc(1,1,inpp), imt, jmt)
call setbcx (sbc(1,1,isr), imt, jmt)
# endif
# else
sbc(:,:,iro) = sbc(:,:,iro)*f1a
# endif
call setbcx (sbc(1,1,iro), imt, jmt)
# if defined O_embm
!-----------------------------------------------------------------------
! zero diagnostic for river discharge and call river model
!-----------------------------------------------------------------------
disch(:,:) = 0.
call rivmodel
# if defined O_sed && !defined O_sed_uncoupled
!-----------------------------------------------------------------------
! apply CaCO3 weathering flux proportional to discharge
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------
# if defined O_sed_constrain_weath
if (weathflx .lt. 0.) weathflx = 0.
# endif
# if defined O_save_carbon_totals
dicwflx = weathflx
# endif
tmp = 0.
do j=2,jmtm1
do i=2,imtm1
if (kmt(i,j) .ne. 0.) then
tmp = tmp + disch(i,j)*dxt(i)*dyt(j)*cst(j)
endif
enddo
enddo
tmp = weathflx/tmp
do j=2,jmtm1
# if defined O_global_sums || defined O_save_carbon_totals
fgs = dyt(j)*cst(j)*segtim/secday
# endif
do i=2,imtm1
if (kmt(i,j) .ne. 0.) then
# if defined O_carbon
sbc(i,j,idicflx) = sbc(i,j,idicflx) + disch(i,j)*tmp
# if defined O_npzd_alk
sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + disch(i,j)*2.*tmp
# endif
# if defined O_save_carbon_totals
carblith = carblith - disch(i,j)*dxt(i)*tmp*fgs
& - sbc(i,j,ibdicfx)*dxt(i)*fgs*dzt(kmt(i,j))
# endif
# if defined O_global_sums
dtoic = dtoic - disch(i,j)*dxt(i)*tmp*fgs
# endif
# endif
endif
enddo
enddo
# endif
# endif
# if defined O_sealev_data_transient && defined O_sealev_salinity
!-----------------------------------------------------------------------
! add in flux from sea level change
!-----------------------------------------------------------------------
do j=2,jmtm1
do i=2,imtm1
if (kmt(i,j) .gt. 0)
& sbc(i,j,isflx) = sbc(i,j,isflx) + fs*dsealev
enddo
enddo
# endif
# if defined O_mom && defined O_embm
call coa
# if defined O_fwa
!-----------------------------------------------------------------------
! add additional fresh water flux anomaly
!-----------------------------------------------------------------------
time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
if (time .ge. fwayri .and. time .le. fwayrf) then
if (areafwa .gt. 0) then
! fwaflxi is in Sv (1e6 m3 s-1) and fwarate is in Sv/yr
fwaflx = fwaflxi + (time - fwayri)*fwarate
! convert to flux in g salt cm-2 s-1
fwaflx = -socn*fwaflx*1.e12/areafwa
# if defined O_fwa_precip
call areaavg (precip, fwawt, tmp)
if (tmp .gt. 0) fwaflx = fwaflx/tmp
# endif
do j=2,jmtm1
do i=2,imtm1
sbc(i,j,isflx) = sbc(i,j,isflx) + fwaflx*fwawt(i,j)
# if defined O_fwa_precip
& *precip(i,j)
# endif
enddo
enddo
endif
if (compensate .and. areafwc .gt. 0) then
! fwaflxi is in Sv (1e6 m3 s-1) and fwarate is in Sv/yr
fwcflx = fwaflxi + (time - fwayri)*fwarate
! convert to opposite flux in g salt cm-2 s-1
fwcflx = socn*fwcflx*1.e12/areafwc
# if defined O_fwa_compevap
call areaavg (evap, fwcwt, tmp)
if (tmp .gt. 0) fwcflx = fwcflx/tmp
# endif
do j=2,jmtm1
do i=2,imtm1
sbc(i,j,isflx) = sbc(i,j,isflx) + fwcflx*fwcwt(i,j)
# if defined O_fwa_compevap
& *evap(i,j)
# endif
enddo
enddo
endif
endif
# endif
# if defined O_carbon || defined O_npzd_alk || defined O_npzd_o2 || defined O_npzd || defined O_cfcs_data || defined O_cfcs_data_transient
!-----------------------------------------------------------------------
! add normalized virtual fluxes to other tracers
!-----------------------------------------------------------------------
tarea = 0.
tsflx = 0.
rsocn = 1./socn
do j=2,jmtm1
do i=2,imtm1
if (tmsk(i,j) .ge. 0.5) then
area = dxt(i)*dyt(j)*cst(j)
tarea = tarea + area
tsflx = tsflx + sbc(i,j,isflx)*area
endif
enddo
enddo
tsflx = tsflx/tarea
do j=2,jmtm1
do i=2,imtm1
if (tmsk(i,j) .ge. 0.5) then
tmp = (sbc(i,j,isflx) - tsflx)*rsocn
vflux(i,j) = tmp
# if defined O_carbon
sbc(i,j,idicflx) = sbc(i,j,idicflx) + gaost(idic)*tmp
# if defined O_carbon_14
sbc(i,j,ic14flx) = sbc(i,j,ic14flx) + gaost(ic14)*tmp
# endif
# endif
# if defined O_npzd_alk
sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + gaost(ialk)*tmp
# endif
# if defined O_npzd_o2
sbc(i,j,io2flx) = sbc(i,j,io2flx) + gaost(io2)*tmp
# endif
# if defined O_npzd
sbc(i,j,ipo4flx) = sbc(i,j,ipo4flx) + gaost(ipo4)*tmp
# if !defined O_npzd_no_vflux
sbc(i,j,iphytflx) = sbc(i,j,iphytflx) + gaost(iphyt)*tmp
sbc(i,j,izoopflx) = sbc(i,j,izoopflx) + gaost(izoop)*tmp
sbc(i,j,idetrflx) = sbc(i,j,idetrflx) + gaost(idetr)*tmp
# endif
# if defined O_npzd_nitrogen
sbc(i,j,ino3flx) = sbc(i,j,ino3flx) + gaost(ino3)*tmp
# if !defined O_npzd_no_vflux
sbc(i,j,idiazflx) = sbc(i,j,idiazflx) + gaost(idiaz)*tmp
# endif
# endif
# endif
# if defined O_cfcs_data || defined O_cfcs_data_transient
sbc(i,j,icfc11flx) = sbc(i,j,icfc11flx) + gaost(icfc11)*tmp
sbc(i,j,icfc12flx) = sbc(i,j,icfc12flx) + gaost(icfc12)*tmp
# endif
endif
enddo
enddo
# if defined O_npzd_alk
# endif
# if defined O_carbon
call setbcx (sbc(1,1,idicflx), imt, jmt)
# if defined O_carbon_14
call setbcx (sbc(1,1,ic14flx), imt, jmt)
# endif
# endif
# if defined O_npzd_alk
call setbcx (sbc(1,1,ialkflx), imt, jmt)
# endif
# if defined O_npzd_o2
call setbcx (sbc(1,1,io2flx), imt, jmt)
# endif
# if defined O_npzd
call setbcx (sbc(1,1,ipo4flx), imt, jmt)
# if !defined O_npzd_no_vflux
call setbcx (sbc(1,1,iphytflx), imt, jmt)
call setbcx (sbc(1,1,izoopflx), imt, jmt)
call setbcx (sbc(1,1,idetrflx), imt, jmt)
# endif
# if defined O_npzd_nitrogen
call setbcx (sbc(1,1,ino3flx), imt, jmt)
# if !defined O_npzd_no_vflux
call setbcx (sbc(1,1,idiazflx), imt, jmt)
# endif
# endif
# endif
# if defined O_cfcs_data || defined O_cfcs_data_transient
call setbcx (sbc(1,1,icfc11flx), imt, jmt)
call setbcx (sbc(1,1,icfc12flx), imt, jmt)
# endif
# endif
# endif
# if defined O_embm_read_sflx
call read_var ("F_salt.nc", "F_salt", sbc(1,1,isflx), ntrec)
# elif defined O_embm_write_sflx
call write_var ("F_salt.nc", "F_salt", sbc(1,1,isflx), ntrec)
# endif
#endif
# if defined O_embm && defined O_sulphate_data || defined O_sulphate_data_transient
!-----------------------------------------------------------------------
! set anthropogenic sulphate forcing (increase in surface albedo)
! for the next atmospheric step from updated surface albedo
!-----------------------------------------------------------------------
call sulphdata
zero(1) = 0.
cz(1) = 0.
! full direct and indirect sulphate forcing
sulphfac = 0.29
# if defined O_sulphate_data_direct && !defined O_sulphate_data_indirect
! reduce sulphate forcing to just the direct effect (0.4/1.1)
sulphfac = sulphfac*.4/1.1
# endif
# if defined O_sulphate_data_indirect && !defined O_sulphate_data_direct
! reduce sulphate forcing to just the indirect effect (0.7/1.1)
sulphfac = sulphfac*.7/1.1
# endif
do j=jsp1,jem1
do i=isp1,iem1
if (sulph(i,j,2) .gt. 0. .and. solins(i,j) .gt. 0.) then
tot = 0.
twt = 0.
! calculate sunlight weighted daily average zenith angle
do k=1,24
call zenith (1, (k-1)*3600., 3600., daylen, tlat(i,j)
&, zero(1), sindec, cz(1))
tot = tot + cz(1)*cz(1)
twt = twt + cz(1)
enddo
cz(1) = 0.
if (twt .gt. 0) cz(1) = tot/twt
! convert from optical depth to albedo
if (cz(1) .gt. 0.05) then
! reonstruct surface coalbedo
tmp = solins(i,j) - solins(i,j)*(1. - sbc(i,j,iaca))
& - solins(i,j)*sbc(i,j,iaca)*scatter
if (tmp .gt. 0.) then
tmp = dnswr(i,j)/tmp
sulph(i,j,2) = sulphfac*sulph(i,j,2)*tmp**2/cz(1)
else
sulph(i,j,2) = 0.
endif
else
sulph(i,j,2) = 0.
endif
endif
enddo
enddo
# endif
return
end
#if defined O_embm_read_sflx
subroutine read_var (fname, vname, var, ntrec)
implicit none
character (*) :: fname, vname
integer ib(10), ic(10), iou, ntrec
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "tmngr.h"
real var(imt,jmt), tmpij(imtm2,jmtm2)
ntrec = ntrec + 1
call openfile (trim(fname), iou)
ib(:) = 1
ib(3) = ntrec
ic(:) = 1
ic(1) = imtm2
ic(2) = jmtm2
call getvara (trim(vname), iou, imtm2*jmtm2, ib, ic, tmpij
&, 1., 0.)
var(2:imtm1,2:jmtm1) = tmpij(1:imtm2:1:jmtm2)
embmbc (var)
return
end
#endif
#if defined O_embm_write_sflx
subroutine write_var (fname, vname, var, ntrec)
implicit none
character (*) :: fname, vname
integer ib(10), ic(10), it(10), id_time, id_xt, id_yt, iou, ntrec
integer nyear, nmonth, nday, nhour, nmin, nsec
real tmp
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "coord.h"
include "tmngr.h"
real var(imt,jmt), tmpij(imtm2,jmtm2), tmpi(imtm2), tmpj(jmtm2)
ntrec = ntrec + 1
call openfile (trim(fname), iou)
if (ntrec .eq. 1) then
call redef (iou)
call defdim ('time', iou, 0, id_time)
call defdim ('longitude', iou, imtm2, id_xt)
call defdim ('latitude', iou, jmtm2, id_yt)
it(:) = id_time
call defvar ('time', iou, 1, it, c0, c0, 'T', 'D'
&, 'time', 'time', 'years since 0-1-1')
call putatttext (iou, 'time', 'calendar', calendar)
it(1) = id_xt
call defvar ('longitude', iou, 1, it, 0., 0., 'X', 'D'
&, 'longitude', 'longitude', 'degrees_east')
it(1) = id_yt
call defvar ('latitude', iou, 1, it, 0., 0., 'Y', 'D'
&, 'latitude', 'latitude', 'degrees_north')
it(:) = id_xt
it(2) = id_yt
it(3) = id_time
call defvar (trim(vname), iou, 3, it, -1.e20, 1.e20, ' ', 'D'
&, ' ',' ', ' ')
call enddef (iou)
!-----------------------------------------------------------------------
! write 1d data (t)
!-----------------------------------------------------------------------
call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
ib(:) = 1
ic(:) = imtm2
tmpi(1:imtm2) = xt(2:imtm1)
call putvara ('longitude', iou, imtm2, ib, ic, tmpi, c1, c0)
ic(:) = jmtm2
tmpj(1:jmtm2) = xt(2:jmtm1)
call putvara ('latitude', iou, jmtm2, ib, ic, tmpj, c1, c0)
endif
call putvars ('time', iou, ntrec, relyr, c1, c0)
ib(:) = 1
ib(3) = ntrec
ic(:) = 1
ic(1) = imtm2
ic(2) = jmtm2
tmpij(1:imtm2:1:jmtm2) = var(2:imtm1,2:jmtm1)
call putvara (trim(vname), iou, imtm2*jmtm2, ib, ic, tmpij
&, 1., 0.)
return
end
#endif