!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module sw_absorption 4,15
!BOP
! !MODULE: sw_absorption
!
! !DESCRIPTION:
! This module computes shortwave absorption
!
! sw_absorption_type
! = 'top-layer' : all sw absorbed in top ocean layer
! = 'jerlov' : sw absorption by jerlov water type
! = 'chlorophyll' : sw absorption by chlorophyll amount
!
! For chlorophyll based absorption, transmission from surface
! (z=0) to level z = A_1 exp(-B_1 z) + A_2 exp(-B_2 z)
!
! A,B coefficients from Table 1a, J. Carter Ohlmann, 2002, Ocean Radiant
! Heating in Climate Models, Journal of Climate, 2003, in press.
! Chlorophyll concentration in mg/m^3. Values for chl < .01 mg/m^3 are
! linear extrpolations from two smallest chlorophyll amounts. Values
! over 3 mg/m^3 up to 10 mg/m^3 were provided by Carter specifically
! for this application. Chlorophyll amounts limited from .001 to 10 mg/m^3
!
! Note that A,B values are relative to the net shortwave flux just below
! the ocean surface [E_d(0-) in Carter notation]. This means ocean
! albedo effects have been accounted for (i.e. E_d(0-) = E_d(0+)(1-alpha_o)
! where E_d(0+) is the down incident shortwave flux at the ocean surface
! and alpha_o is the ocean albedo).
!
! Note also that the second order effects of solar zenith angle and
! cloud cover have been ignored in this implementation.
!
! Chlorophyll monthly distribution is provided over all ocean, and used to
! compute regional distribution of shortwave absorption in the upper ocean.
!
! Chlorophyll based transmissions are calculated once a month and
! held fixed over the month in this implementation.
! !REVISION HISTORY:
! SVN:$Id: sw_absorption.F90 22832 2010-05-07 23:05:11Z njn01 $
! !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_tools
use time_management
use prognostic
use forcing_shf
use tavg, only: define_tavg_field, tavg_requested, accumulate_tavg_field
use exit_mod
use registry, only: registry_match
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_sw_absorption, &
add_sw_absorb, &
sw_absorb_frac, &
set_chl, &
sw_trans_chl
! !PUBLIC DATA MEMBERS:
character (char_len), public :: &
sw_absorption_type ! type of sw_absorption
real (r8), dimension(nx_block,ny_block,max_blocks_clinic),public :: &
TRANS, &! transmission (1 to 0) from surface
! to depth ZTRANS
TRANSKM1, &! transmission from surface to level k-1
ZTRANS ! depth of transmission cm
!EOP
!BOC
character (char_len) :: &
chl_option, &! chlorophyll option ['file', 'model']
chl_filename, &! chlorophyll data file name
chl_file_fmt, &! chlorophyll data file format
chl_data_name ! chlorophyll short name (eg, 'CHL')
integer (int_kind) :: &
chl_bndy_loc, &! location and field types for ghost
chl_bndy_type ! cell update routines
integer (int_kind) :: &
jerlov_water_type ! jerlov water type from 1 to 5
real (r8), dimension(0:km) :: &
sw_absorb ! sw transmission for jerlov water type
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
CHL ! chlorophyll amount mg/m^3
real (r8) :: &
chlmin, &! minimum chlorophyll amount mg/m^3
chlmax, &! maximum chlorophyll amount mg/m^3
dlogchl ! delta log chl
integer (int_kind), dimension(nx_block,ny_block,max_blocks_clinic) :: &
CHLINDX ! chlorophyll index from table
real (r8), allocatable, dimension(:,:,:,:) :: &
CHL_DATA ! global chlorophyll climatology by month
integer (int_kind), parameter :: &
nsub = 400, &! total number chlorophyll concentrations in look up table
ksol = 2*km, &! (0:ksol) transmission levels (layer interface + mid-point)
nchl = 31 ! number of chlorophyll absorption coefficients
real (r8), dimension(0:ksol,0:nsub) :: &
Tr ! table look-up transmissions (level x chlorophyll)
real (r8), dimension(0:ksol) :: &
ztr ! transmission array levels (layer interface + mid-point)
real (r8) :: &
A1,B1,A2,B2 ! A,B coefficient interpolation values
real (r8), parameter, dimension(nchl) :: &!chl for table look-up
chlcnc = (/ &
.001_r8, .005_r8, .01_r8, .02_r8, &
.03_r8, .05_r8, .10_r8, .15_r8, &
.20_r8, .25_r8, .30_r8, .35_r8, &
.40_r8, .45_r8, .50_r8, .60_r8, &
.70_r8, .80_r8, .90_r8, 1.00_r8, &
1.50_r8, 2.00_r8, 2.50_r8, 3.00_r8, &
4.00_r8, 5.00_r8, 6.00_r8, 7.00_r8, &
8.00_r8, 9.00_r8,10.00_r8 /)
real (r8), parameter, dimension(nchl) :: &
A_1 = (/ &
0.4421_r8, 0.4451_r8, 0.4488_r8, &
0.4563_r8, &
0.4622_r8, 0.4715_r8, 0.4877_r8, &
0.4993_r8, &
0.5084_r8, 0.5159_r8, 0.5223_r8, &
0.5278_r8, &
0.5326_r8, 0.5369_r8, 0.5408_r8, &
0.5474_r8, &
0.5529_r8, 0.5576_r8, 0.5615_r8, &
0.5649_r8, &
0.5757_r8, 0.5802_r8, 0.5808_r8, &
0.5788_r8, &
0.56965_r8,0.55638_r8,0.54091_r8, &
0.52442_r8, &
0.50766_r8,0.49110_r8,0.47505_r8 /)
real (r8), parameter, dimension(nchl) :: &
A_2 = (/ &
0.2981_r8, 0.2963_r8, 0.2940_r8, &
0.2894_r8, &
0.2858_r8, 0.2800_r8, 0.2703_r8, &
0.2628_r8, &
0.2571_r8, 0.2523_r8, 0.2481_r8, &
0.2444_r8, &
0.2411_r8, 0.2382_r8, 0.2356_r8, &
0.2309_r8, &
0.2269_r8, 0.2235_r8, 0.2206_r8, &
0.2181_r8, &
0.2106_r8, 0.2089_r8, 0.2113_r8, &
0.2167_r8, &
0.23357_r8,0.25504_r8,0.27829_r8, &
0.30274_r8, &
0.32698_r8,0.35056_r8,0.37303_r8 /)
real (r8), parameter, dimension(nchl) :: &
B_1 = (/ &
0.0287_r8, 0.0301_r8, 0.0319_r8, &
0.0355_r8, &
0.0384_r8, 0.0434_r8, 0.0532_r8, &
0.0612_r8, &
0.0681_r8, 0.0743_r8, 0.0800_r8, &
0.0853_r8, &
0.0902_r8, 0.0949_r8, 0.0993_r8, &
0.1077_r8, &
0.1154_r8, 0.1227_r8, 0.1294_r8, &
0.1359_r8, &
0.1640_r8, 0.1876_r8, 0.2082_r8, &
0.2264_r8, &
0.25808_r8,0.28498_r8,0.30844_r8, &
0.32932_r8, &
0.34817_r8,0.36540_r8,0.38132_r8 /)
real (r8), parameter, dimension(nchl) :: &
B_2 = (/ &
0.3192_r8, 0.3243_r8, 0.3306_r8, &
0.3433_r8, &
0.3537_r8, 0.3705_r8, 0.4031_r8, &
0.4262_r8, &
0.4456_r8, 0.4621_r8, 0.4763_r8, &
0.4889_r8, &
0.4999_r8, 0.5100_r8, 0.5191_r8, &
0.5347_r8, &
0.5477_r8, 0.5588_r8, 0.5682_r8, &
0.5764_r8, &
0.6042_r8, 0.6206_r8, 0.6324_r8, &
0.6425_r8, &
0.66172_r8,0.68144_r8,0.70086_r8, &
0.72144_r8, &
0.74178_r8,0.76190_r8,0.78155_r8 /)
integer (int_kind) :: &
tavg_QSW_HTP, & ! tavg id for QSW_HTP (solar short-wave heat flux in top layer)
tavg_QSW_3D ! tavg id for 3D QSW at top of cell
!-----------------------------------------------------------------------
! named field indices
!-----------------------------------------------------------------------
integer (kind=int_kind) :: chl_nf_ind = 0
!EOC
!***********************************************************************
contains
!***********************************************************************
subroutine init_sw_absorption 1,27
!-----------------------------------------------------------------------
!
! initialize shortwave absorption.
!
!-----------------------------------------------------------------------
use named_field_mod, only: named_field_get_index
implicit none
integer (int_kind) :: &
k,n, &! level index
nu, &! unit for input dataset
nml_error, &! namelist error flag
bid, &! local block address for this block
iblock
real (r8), allocatable, dimension(:,:,:,:) :: &
TEMP_DATA ! temporary data array
type (datafile) :: &
forcing_file !data file structure for input forcing file
type (io_field_desc) :: &
io_chl ! io field descriptor for input sss field
type (io_dim) :: &
i_dim, j_dim, &! dimension descriptors for horiz dimensions
month_dim ! dimension descriptor for monthly data
namelist /sw_absorption_nml/ &
sw_absorption_type, &! 'top-layer', 'jerlov' or 'chlorophyll'
jerlov_water_type, &! jerlov water type from 1 to 5
chl_option, &! chlorophyll option ['file', 'model']
chl_filename, &! include local filepath in name
chl_file_fmt ! include local filepath in name
!-----------------------------------------------------------------------
if (.not. registry_match('init_passive_tracers')) then
call exit_POP(sigAbort, 'init_passive_tracers not called ' /&
&/ 'before init_sw_absorption. This is necessary to ' /&
&/ 'allow for registry of model_chlorophyll.')
end if
!-----------------------------------------------------------------------
!
! read shortwave absorption namelist input after setting default values.
!
!-----------------------------------------------------------------------
sw_absorption_type = 'jerlov'
jerlov_water_type = 3
chl_option = 'file'
chl_filename = 'unknown-chl'
chl_file_fmt = 'bin'
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=sw_absorption_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 sw_absorption_nml')
endif
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,blank_fmt)
write(stdout,*) ' Short-wave absorption:'
write(stdout,blank_fmt)
write(stdout,*) ' sw_absorption namelist '
write(stdout,blank_fmt)
write(stdout, sw_absorption_nml)
write(stdout,blank_fmt)
endif
call broadcast_scalar(sw_absorption_type, master_task)
call broadcast_scalar(jerlov_water_type, master_task)
call broadcast_scalar(chl_option, master_task)
call broadcast_scalar(chl_filename, master_task)
if (sw_absorption_type .ne. 'top-layer'.and. &
sw_absorption_type .ne. 'jerlov'.and. &
sw_absorption_type .ne. 'chlorophyll' ) then
call exit_POP(sigAbort,'ERROR sw_absorption_type unknown')
endif
if (sw_absorption_type .eq. 'chlorophyll') then
if (chl_option .ne. 'file' .and. chl_option .ne. 'model') then
call exit_POP(sigAbort,'ERROR chl_option unknown')
endif
if (chl_option .eq. 'model') then
call named_field_get_index('model_chlorophyll', chl_nf_ind, &
exit_on_err=.false.)
if (chl_nf_ind == 0) then
call exit_POP(sigAbort,'chl_option==model, but ' /&
&/ 'model_chlorophyll is not registered')
endif
endif
endif
!-----------------------------------------------------------------------
! define shortwave solar absorption model.
!-----------------------------------------------------------------------
select case (sw_absorption_type)
case ('top-layer')
sw_absorb(0) = c1
do k = 1, km
sw_absorb(k) = c0
enddo
case ('jerlov')
sw_absorb(0) = c1
sw_absorb(km) = c0
do k = 1, km-1
call sw_absorb_frac(zw(k),sw_absorb(k))
enddo
case ('chlorophyll')
call set_chl_trn ! compute transmission table
if (chl_option == 'file') then
!-----------------------------------------------------------------------
!
! Read in chlorophyll data
!
!-----------------------------------------------------------------------
chl_bndy_loc = field_loc_center
chl_bndy_type = field_type_scalar
chl_data_name = 'CHL'
allocate( CHL_DATA (nx_block,ny_block,max_blocks_clinic,12))
allocate( TEMP_DATA(nx_block,ny_block,12,max_blocks_clinic))
CHL_DATA = c0
forcing_file = construct_file(chl_file_fmt, &
full_name=trim(chl_filename), &
record_length = rec_type_dbl, &
recl_words=nx_global*ny_global)
call data_set(forcing_file,'open_read')
i_dim = construct_io_dim('i',nx_global)
j_dim = construct_io_dim('j',ny_global)
month_dim = construct_io_dim('month',12)
io_chl = construct_io_field( &
trim(chl_data_name), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = chl_bndy_loc, &
field_type = chl_bndy_type, &
d3d_array=TEMP_DATA(:,:,:,:))
call data_set(forcing_file,'define',io_chl)
call data_set(forcing_file,'read' ,io_chl)
call destroy_io_field(io_chl)
!*** re-order data
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
CHL_DATA(:,:,iblock,n) = TEMP_DATA(:,:,n,iblock)
end do
end do
!$OMP END PARALLEL DO
deallocate(TEMP_DATA)
call data_set(forcing_file,'close')
call destroy_file(forcing_file)
if (my_task.eq.master_task) then
write(stdout,blank_fmt)
write(stdout,*) ' CHL monthly file read: ',chl_filename
endif
!-----------------------------------------------------------------------
!
! Set initial chlorophyll amount
!
!-----------------------------------------------------------------------
CHL(:,:,:) = CHL_DATA(:,:,:,imonth) ! constant across month
CHL = max(CHL,chlmin) ! set lowest allowed limit
CHL = min(CHL,chlmax) ! set highest allowed limit
CHLINDX = log10(CHL/chlmin)/dlogchl ! compute chlorphyll index
CHLINDX = max(CHLINDX,0) ! minimum limit for chl index
CHLINDX = min(CHLINDX,nsub) ! maximum limit for chl index
end if ! chl_option == 'file'
end select
call define_tavg_field(tavg_QSW_HTP,'QSW_HTP',2, &
long_name='Solar Short-Wave Heat Flux in top layer', &
units='watt/m^2', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_QSW_3D,'QSW_3D',3, &
long_name='Solar Short-Wave Heat Flux', &
units='watt/m^2', grid_loc='3112', &
coordinates='TLONG TLAT z_w_top time')
!-----------------------------------------------------------------------
end subroutine init_sw_absorption
!***********************************************************************
subroutine set_chl 1,3
use named_field_mod, only: named_field_get
logical (kind=log_kind) :: first_call = .true.
integer (int_kind) :: &
iblock
!-----------------------------------------------------------------------
!
! check that set_sflux_passive_tracers has been called
!
!-----------------------------------------------------------------------
if (first_call .and. chl_option == 'model') then
if (.not. registry_match('set_sflux_passive_tracers')) then
call exit_POP(sigAbort, 'set_sflux_passive_tracers not ' /&
&/ 'called before set_chl. This is necessary for ' /&
&/ 'exact restart when Chl is prognostic.')
end if
endif
!-----------------------------------------------------------------------
!
! Update chlorophyll amount; update every new month
!
!-----------------------------------------------------------------------
if( sw_absorption_type .eq. 'chlorophyll' ) then
if (imonth .ne. imonth_last .or. & ! if new month, update chl
chl_option == 'model' ) then ! update every timestep for 'model'
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
if (chl_option == 'file') then
CHL(:,:,iblock) = CHL_DATA(:,:,iblock,imonth) ! constant across month
else
call named_field_get(chl_nf_ind, iblock, CHL(:,:,iblock))
endif
CHL(:,:,iblock) = max(CHL(:,:,iblock),chlmin) ! set lowest allowed limit
CHL(:,:,iblock) = min(CHL(:,:,iblock),chlmax) ! set highest allowed limit
CHLINDX(:,:,iblock) = log10(CHL(:,:,iblock)/chlmin)/dlogchl ! compute chlorphyll index
CHLINDX(:,:,iblock) = max(CHLINDX(:,:,iblock),0) ! minimum limit for chl index
CHLINDX(:,:,iblock) = min(CHLINDX(:,:,iblock),nsub) ! maximum limit for chl index
end do
!$OMP END PARALLEL DO
endif
endif
first_call = .false.
end subroutine set_chl
!***********************************************************************
subroutine set_chl_trn 1,3
!-----------------------------------------------------------------------
!
! Set chlorophyll transmissions for table look-up.
!
! Look-up table coefficients computed over the range .001
! to 10 mg/m^3, by a constant step in log(chl). We
! choose to have 100 steps per decade, so increment
! in chl is around 2.3% .
!
! Note that table coefficients require depth in m, so cm to m
! unit conversion is done where appropriate.
!
!-----------------------------------------------------------------------
implicit none
real (r8) :: &
chlamnt, &! chlorophyll amount (density) in mg/m^3
w1, &! interpolation weight
w2 ! interpolation weight
real (r8) :: &
logchl, &! log of chlorophyll
zprnt ! z level for print
integer (int_kind) :: &
n, &! index for absorber amount
m, &! index for absorber amount
mc, &! specific absorber index
k, &! vertical index
kint, &! vertical index for layer interface
kmid ! vertical index for layer mid-point
real (r8) :: &
Trprnt, &! chlorophyll transmissions for diagnostic print
Trjrlv, &! Jerlov transmissions for diagnostic print
arg ! argument for exponential
real (r8), parameter :: &
maxarg = 35.0_r8 ! maximum allowed exponential argument
logical (log_kind) :: &
prnt = .false. ! if true, generates diagnostic prints
!-----------------------------------------------------------------------
!
! Diagnostic prints if desired
!
!-----------------------------------------------------------------------
if (my_task == master_task) then
if (prnt) then
write(stdout,5)
5 format( &
' Transmissions for selected chlorphyll amounts'/ &
' Interface Mid-point Transmission Jerlov Tr (type 3 = IB)')
do n=1,nchl
if (n .eq. 1 .or. n .eq. 3 .or. &
n .eq. 6 .or. n .eq. 7 .or. &
n .eq. 15 .or. n .eq. 20 .or. &
n .eq. 24 .or. n .eq. 31) then
write(stdout,8) chlcnc(n)
8 format(/' chlorophyll amount = ',f6.3/)
do k=1,km
kint = k-1
if (kint.eq.0 ) then
Trprnt = c1
else
arg = min(B_1(n)*zw(kint)*mpercm,maxarg)
Trprnt = A_1(n)*exp(-arg)
arg = min(B_2(n)*zw(kint)*mpercm,maxarg)
Trprnt = Trprnt + A_2(n)*exp(-arg)
endif
arg = min(zw(kint)*mpercm/1.0_r8,maxarg)
Trjrlv = .67_r8*exp(-arg)
arg = min(zw(kint)*mpercm/17.0_r8,maxarg)
Trjrlv = Trjrlv + .33_r8*exp(-arg)
if (kint .eq. 0 ) then
zprnt = c0
Trjrlv = c1
if (Trprnt .gt. .00005 ) &
write(stdout,10) kint,zprnt,Trprnt,Trjrlv
10 format(1x,i3,1x,f6.1,15x,f6.4,8x,f6.4)
else
if (Trprnt .gt. .00005 ) &
write(stdout,11) kint,zw(kint)*mpercm,Trprnt,Trjrlv
11 format(1x,i3,1x,f6.1,16x,f5.4,8x,f5.4)
endif
if (k .le. km-1 ) then
kmid = kint + 1
arg = min(B_1(n)*zt(kmid)*mpercm,maxarg)
Trprnt = A_1(n)*exp(-arg)
arg = min(B_2(n)*zt(kmid)*mpercm,maxarg)
Trprnt = Trprnt + A_2(n)*exp(-arg)
arg = min(zt(kmid)*mpercm/1.0_r8,maxarg)
Trjrlv = .67_r8*exp(-arg)
arg = min(zt(kmid)*mpercm/17.0_r8,maxarg)
Trjrlv = Trjrlv + .33_r8*exp(-arg)
if (Trprnt .gt. .00005 ) &
write(stdout,15) kmid,zt(kmid)*mpercm,Trprnt,Trjrlv
15 format(12x,i3,1x,f6.1,5x,f5.4,8x,f5.4)
endif
enddo
endif
enddo
endif
endif
!-----------------------------------------------------------------------
!
! Compute transmission array
!
!-----------------------------------------------------------------------
ztr(0) = c0
do k=1,km
ztr(2*k-1) = zt(k)
ztr(2*k) = zw(k)
enddo
chlmin = chlcnc(1)
chlmax = chlcnc(nchl)
dlogchl = (log10(chlmax)-log10(chlmin))/real(nsub)
if (my_task == master_task) then
if (prnt) then
write(stdout,20) chlmin,chlmax,nsub,dlogchl
20 format(/1x,' minimum chlorophyll = ',1pe11.4,' mg/m^3'/ &
1x,' maximum chlorophyll = ',1pe11.4,' mg/m^3'/ &
1x,' number of intervals = ',i5/ &
1x,' log chlorophyll step = ',1pe11.4)
endif
endif
logchl = log10(chlmin) - dlogchl
do n=0,nsub
logchl = logchl + dlogchl
chlamnt = 10**(logchl)
do m=1,nchl-1
if (chlcnc(m) .le. chlamnt .and. &
chlamnt .le. chlcnc(m+1) ) then
mc = m
goto 30
endif
enddo
if (my_task == master_task) then
write(stdout,25) chlamnt
25 format(/' Could not find range for chlamnt = ',1pe11.4)
call exit_POP(sigAbort,' set_chl_trans range error for chlamnt ')
endif
30 continue
w2 = (chlamnt-chlcnc(mc))/(chlcnc(mc+1)-chlcnc(mc))
w1 = c1 - w2
A1 = A_1(mc)*w1 + A_1(mc+1)*w2
A2 = A_2(mc)*w1 + A_2(mc+1)*w2
B1 = B_1(mc)*w1 + B_1(mc+1)*w2
B2 = B_2(mc)*w1 + B_2(mc+1)*w2
if (my_task == master_task) then
if (prnt) then
write(stdout,35) n,chlamnt,A1,A2,B1,B2
35 format(/' n, chl = ',i5,2x,1pe11.4/ &
' A1,A2,B1,B2 = ',4(f6.4,2x)/ &
' Level Depth(m) Transmission ')
endif
endif
Tr(0,n) = c1
do k=1,ksol
arg = min(B1*ztr(k)*mpercm,maxarg)
Tr(k,n) = A1*exp(-arg)
arg = min(B2*ztr(k)*mpercm,maxarg)
Tr(k,n) = Tr(k,n) + A2*exp(-arg)
enddo
if (my_task == master_task) then
if (prnt) then
do k=0,ksol
if (Tr(k,n) .gt. .00005 ) &
write(stdout,40) k,ztr(k)*mpercm,Tr(k,n)
40 format(1x,i3,4x,f6.1,5x,f6.4)
enddo
endif
endif
enddo
if (my_task.eq.master_task) then
write(stdout,blank_fmt)
write(stdout,*) ' Chlorophyll transmission table computed'
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
end subroutine set_chl_trn
!***********************************************************************
!BOP
! !IROUTINE: sw_absorb_frac
! !INTERFACE:
subroutine sw_absorb_frac( depth, sw_absorb_fraction ) 4
! !DESCRIPTION:
! Computes fraction of solar short-wave flux penetrating to
! specified depth due to exponential decay in Jerlov water type.
! Reference : two band solar absorption model of Simpson and
! Paulson (1977)
! Note: below 200m the solar penetration gets set to zero,
! otherwise the limit for the exponent ($+/- 5678$) needs to be
! taken care of.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8) :: &
depth ! vertical depth (cm, >0.) for desired sw fraction
! !OUTPUT PARAMETERS:
real (r8) :: &
sw_absorb_fraction ! short wave (radiation) fractional decay
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
num_water_types = 5 ! max number of different water types
real (r8), parameter :: &
depth_cutoff = -200.0_r8
real (r8) :: &
depth_neg_meters
!-----------------------------------------------------------------------
!
! define Jerlov water properties with rfac, depth1, depth2
! Jerlov water type : I IA IB II III
! jerlov_water_type : 1 2 3 4 5
!
!-----------------------------------------------------------------------
real (r8), dimension(num_water_types) :: &
rfac = (/ 0.58_r8, 0.62_r8, 0.67_r8, 0.77_r8, 0.78_r8 /), &
depth1 = (/ 0.35_r8, 0.60_r8, 1.00_r8, 1.50_r8, 1.40_r8 /), &
depth2 = (/ 23.0_r8, 20.0_r8, 17.0_r8, 14.0_r8, 7.90_r8 /)
!-----------------------------------------------------------------------
!
! compute absorption fraction
!
!-----------------------------------------------------------------------
depth_neg_meters = -depth*mpercm ! convert from cm to m and
! change sign
if (depth_neg_meters < depth_cutoff) then
sw_absorb_fraction = c0
else
sw_absorb_fraction = rfac(jerlov_water_type)* &
exp(depth_neg_meters/depth1(jerlov_water_type)) + &
(c1 - rfac(jerlov_water_type))* &
exp(depth_neg_meters/depth2(jerlov_water_type))
endif
!-----------------------------------------------------------------------
!EOC
end subroutine sw_absorb_frac
!***********************************************************************
!BOP
! !IROUTINE: add_sw_absorb
! !INTERFACE:
subroutine add_sw_absorb(T_SOURCE, SHF_QSW, k, this_block) 1,9
! !DESCRIPTION:
! If surface short wave heat flux is available, this routine caculates
! the flux which passes through the top layer and enters lower vertical
! depths in the ocean. This flux is added as a source term in the
! baroclinic equations.
!
! !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)
! sw absorption added only to other potential
! temperature tracers
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k ! vertical level index
type (block), intent(in) :: &
this_block ! block info for this block
real (r8), dimension(nx_block,ny_block), intent(in) :: &
SHF_QSW
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
bid ! local block index
real (r8), dimension(nx_block,ny_block) :: &
WORK ! temporary work space
!-----------------------------------------------------------------------
!
! calculate short wave absorption if available
! absorption profile pre-calculated in init_shf routine
!
!-----------------------------------------------------------------------
if (lsw_absorb .or. registry_match('lcoupled')) then ! short wave flux is available
bid = this_block%local_id
WORK = max(SHF_QSW,c0) !*** insure no neg QSW - store in work
select case (sw_absorption_type)
case ('top-layer', 'jerlov')
if (partial_bottom_cells) then
where (k < KMT(:,:,bid))
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + &
WORK*(sw_absorb(k-1) - sw_absorb(k))* &
dzr(k)
elsewhere ! do not allow energy absorption by the ground
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + WORK* &
(sw_absorb(k-1))/DZT(:,:,k,bid)
endwhere
else
where (k < KMT(:,:,bid))
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + &
WORK*(sw_absorb(k-1) - sw_absorb(k))* &
dzr(k)
elsewhere ! do not allow energy absorption by the ground
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + &
WORK *(sw_absorb(k-1))*dzr(k)
endwhere
endif
if (k == 1) then
if (tavg_requested(tavg_QSW_HTP)) then
call accumulate_tavg_field(WORK*(sw_absorb(0)-sw_absorb(1))/hflux_factor, &
tavg_QSW_HTP,bid,k)
endif
endif
if (tavg_requested(tavg_QSW_3D)) then
call accumulate_tavg_field(WORK*sw_absorb(k-1)/hflux_factor, &
tavg_QSW_3D,bid,k)
endif
case ('chlorophyll')
if( k.eq.1 ) TRANSKM1(:,:,bid) = c1 ! surface z=0
call sw_trans_chl(2*k,this_block) ! use explicit index and
! chl to get transmission
if (partial_bottom_cells) then
where (k < KMT(:,:,bid))
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + &
WORK*(TRANSKM1(:,:,bid)-TRANS(:,:,bid))*dzr(k)
elsewhere ! do not allow energy absorption by the ground
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + &
WORK*TRANSKM1(:,:,bid)/DZT(:,:,k,bid)
endwhere
else
where (k < KMT(:,:,bid))
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + WORK*(TRANSKM1(:,:,bid)- &
TRANS(:,:,bid))*dzr(k)
elsewhere ! do not allow energy absorption by the ground
T_SOURCE(:,:,1) = T_SOURCE(:,:,1) + WORK*TRANSKM1(:,:,bid)*dzr(k)
endwhere
endif
if (k == 1) then
if (tavg_requested(tavg_QSW_HTP)) then
call accumulate_tavg_field(WORK*(TRANSKM1(:,:,bid)-TRANS(:,:,bid))/hflux_factor, &
tavg_QSW_HTP,bid,k)
endif
endif
if (tavg_requested(tavg_QSW_3D)) then
call accumulate_tavg_field(WORK*TRANSKM1(:,:,bid)/hflux_factor, &
tavg_QSW_3D,bid,k)
endif
TRANSKM1(:,:,bid) = TRANS(:,:,bid)
end select
endif
!-----------------------------------------------------------------------
!EOC
end subroutine add_sw_absorb
!***********************************************************************
subroutine sw_trans_chl(kin,this_block) 4
!-----------------------------------------------------------------------
!
! Given level kin (>0) or depth (positive) and chlorophyll amount
! chl, find transmission from table look-up values.
!
!-----------------------------------------------------------------------
implicit none
! !INPUT PARAMETERS:
type (block), intent(in) :: &
this_block ! block info for this block
integer (int_kind), intent(in) :: &
kin ! input index (if > 0), or implied ZTRANS (kin = 0)
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,k, &! horizontal, vertical indices
bid
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
w1, &! weight for level interpolation
w2 ! weight for level interpolation
integer (int_kind), dimension(nx_block,ny_block) :: &
kindx ! interpolation level index
bid = this_block%local_id
!-----------------------------------------------------------------------
!
! Explicit level known
!
!-----------------------------------------------------------------------
if (kin .gt. 0 ) then
do j=1,ny_block
do i=1,nx_block
TRANS(i,j,bid) = Tr(kin ,CHLINDX(i,j,bid))
enddo
enddo
else
!-----------------------------------------------------------------------
!
! Find level indices and interpolate (chlorophyll index
! has previously been computed when chl was updated)
!
!-----------------------------------------------------------------------
kindx = ksol - 1
!-----------------------------------------------------------------------
! F90 compiler requires conditional to be a scalar;
! use explicit indices
!-----------------------------------------------------------------------
do k=1,ksol
do j=1,ny_block
do i=1,nx_block
if (ztr(k-1).le.ZTRANS(i,j,bid) .and. ZTRANS(i,j,bid).lt.ztr(k) ) then
w2(i,j,bid) = (ZTRANS(i,j,bid)-ztr(k-1))/(ztr(k)-ztr(k-1))
w1(i,j,bid) = c1 - w2(i,j,bid)
kindx(i,j) = k - 1
endif
enddo
enddo
enddo
!-----------------------------------------------------------------------
! F90 compiler requires a vector subscript to be a rank 1 expression;
! use explicit indices
!-----------------------------------------------------------------------
do j=1,ny_block
do i=1,nx_block
TRANS(i,j,bid) = w1(i,j,bid)*Tr(kindx(i,j) ,CHLINDX(i,j,bid)) + &
w2(i,j,bid)*Tr(kindx(i,j)+1,CHLINDX(i,j,bid))
enddo
enddo
endif
end subroutine sw_trans_chl
!***********************************************************************
end module sw_absorption
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||