module ms_balance 2,12
!BOP
! !MODULE: ms_balance
! !DESCRIPTION:
! This module contains routines necessary for the balancing of evaporation,
! precipitation, melt, runoff, and salt in marginal-seas regions of the
! ccsm coupled ocean model
!
! !REVISION HISTORY:
! SVN:$Id: ms_balance.F90 21356 2010-03-01 22:12:38Z njn01 $
!
! !USES:
use POP_KindsMod
use POP_IOUnitsMod
use kinds_mod
use domain_size
use domain
use global_reductions
use gather_scatter
use grid
use communicate
use io_tools
use ice
use constants
use time_management
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_ms_balance, &
ms_balancing
! !PUBLIC DATA MEMBERS:
real (r8), dimension(:,:,:,:), allocatable :: MASK_FRAC
real (r8), dimension(max_regions) :: &
annual_depth ,&
monthly_depth
!EOP
!BOC
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_ms_balance
! !INTERFACE:
subroutine init_ms_balance 1,30
! !DESCRIPTION:
! Provides initialization for marginal-sea balancing
!
! For a given marginal sea, determine the corresponding active-ocean
! distribution region for balancing evaporation, precipitation, melt, and runoff
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: max_points = 5000, &
max_iter = 16
logical (log_kind) :: & ! Is this point a:
duplicate ,& ! duplicate point?
marginal_sea ,& ! marginal sea point?
land ,& ! land point?
debug1 = .false. ,& ! Print level-1 debug statements
debug2 = .false. ! Print level-2 debug statements
integer (int_kind) :: &
n ,& ! region number
i,j ,& ! global do-loop indices
npts ,& ! number of redistribution points
iter, nn, pt ,& ! do-loop indices
region_number ,&
mask_index = 0
integer (int_kind), dimension (max_points) :: &
ipts ,& ! storage for i&j points in
jpts ! distribution region
integer (int_kind), allocatable, dimension(:,:) :: &
REGION_MASK_G ! global region mask
real (r8) :: &
mblat, mblon ,& ! center of lat & lon search
lat_reach, lon_reach ,& ! search reach in lat & lon
area ,& ! search area
sum_frac ! sum of distribution fractions
real (r8),dimension (max_points) :: &
areas ,& ! distribution areas
fracs ! distribution fractions
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
WORK
real (r8), allocatable, dimension(:,:) :: &
TLAT_G ,& ! global latitude of cell center
TLON_G ,& ! global longitude of cell center
AREAT_G ,& ! global areas centered at T points
MASK_G ! global redistribution mask
if (my_task == master_task) then
write(stdout,delim_fmt)
write(stdout,blank_fmt)
write(stdout,'(a)') ' Marginal-sea balancing information'
write(stdout,blank_fmt)
write(stdout,delim_fmt)
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
allocate (MASK_FRAC (nx_block,ny_block,max_blocks_clinic,max_ms) )
annual_depth = c0
monthly_depth = c0
!-----------------------------------------------------------------------
! create global arrays of TLON, TLAT, AREAT_G, and REGION_MASK
!-----------------------------------------------------------------------
allocate (REGION_MASK_G(nx_global,ny_global))
allocate (TLON_G (nx_global,ny_global), &
TLAT_G (nx_global,ny_global), &
AREAT_G (nx_global,ny_global), &
MASK_G (nx_global,ny_global) )
WORK = TLON*radian ! TLON in degrees
call gather_global(TLON_G, WORK, master_task,distrb_clinic)
WORK = TLAT*radian ! TLAT in degrees
call gather_global(TLAT_G, WORK, master_task,distrb_clinic)
WORK = DXT*DYT
call gather_global(AREAT_G , WORK, master_task,distrb_clinic)
call gather_global(REGION_MASK_G,REGION_MASK,master_task,distrb_clinic)
if (debug1) call print_regions (REGION_MASK_G)
!-----------------------------------------------------------------------
! for each marginal sea, determine associated distribution region:
!-----------------------------------------------------------------------
marginal_seas: do n=1,num_regions
if (region_info(n)%marginal_sea) then
mask_index = mask_index + 1
MASK_G = c0
!-----------------------------------------------------------------------
! search globally for associated regions
!-----------------------------------------------------------------------
if (my_task == master_task) then
npts = 0
ipts = 0
jpts = 0
areas = c0
fracs = c0
area = c0
mblat = region_info(n)%ms_bal%lat
mblon = region_info(n)%ms_bal%lon
!-----------------------------------------------------------------------
! set initial extent of search in latitude and longitude
!-----------------------------------------------------------------------
lat_reach = c1
lon_reach = c1
iter_loop: do iter = 1, max_iter
!-----------------------------------------------------------------------
! select active-ocean points which lie within search area
!-----------------------------------------------------------------------
j_loop: do j=1,ny_global
i_loop: do i=1,nx_global
if (mblat-lat_reach <= TLAT_G(i,j) .and. &
mblat+lat_reach >= TLAT_G(i,j) .and. &
mblon-lon_reach <= TLON_G(i,j) .and. &
mblon+lon_reach >= TLON_G(i,j) ) then
!-----------------------------------------------------
! is this point identical to a previously selected one?
!-----------------------------------------------------
duplicate = .false.
dup_loop: do nn=1,npts
if (ipts(nn) == i .and. jpts(nn) == j) then
duplicate = .true.
exit dup_loop
endif
end do dup_loop
!-----------------------------------------------------
! is this point a marginal-sea point?
!-----------------------------------------------------
if (REGION_MASK_G(i,j) < 0) then
marginal_sea = .true.
else
marginal_sea = .false.
endif
!-----------------------------------------------------
! is this point a land point?
!-----------------------------------------------------
if (REGION_MASK_G(i,j) == 0) then
land = .true.
else
land = .false.
endif
!--------------------------------------------------
! reject duplicate points
!--------------------------------------------------
if (duplicate) then
if(debug2) call pt_print('reject duplicate point',i,j)
!--------------------------------------------------
! reject marginal-sea points
!--------------------------------------------------
else if (marginal_sea) then
if(debug2) &
call pt_print('reject marginal sea point',i,j)
!--------------------------------------------------
! reject land points
!--------------------------------------------------
else if (land) then
if(debug2) call pt_print('reject land point',i,j)
!--------------------------------------------------
! select unique active-ocean points
!--------------------------------------------------
else
!--------------------------------------------------
! has maximum number of distribution points been
! selected?
!--------------------------------------------------
if (npts == max_points) then
exit iter_loop
endif
npts = npts + 1
ipts (npts) = i
jpts (npts) = j
areas(npts) = AREAT_G(i,j)
area = area + areas(npts)
!---------------------------------------------------
! is search area acceptable size?
!---------------------------------------------------
if( area >= region_info(n)%ms_bal%area)then
if (debug1) &
write(stdout,1002)'(init_ms_balance) ', &
'search area = ', region_info(n)%ms_bal%area
exit iter_loop
endif
endif
endif
end do i_loop
end do j_loop
!------------------------------------------------------------
! increase search area
!------------------------------------------------------------
lat_reach = lat_reach + c1
lon_reach = lon_reach + c1
end do iter_loop
!-----------------------------------------------------------------------
! end of master_task region
!-----------------------------------------------------------------------
endif
!-----------------------------------------------------------------------
! spread the news
!-----------------------------------------------------------------------
call broadcast_scalar(area , master_task)
call broadcast_scalar(npts , master_task)
call broadcast_array (ipts , master_task)
call broadcast_array (jpts , master_task)
call broadcast_array (areas, master_task)
!-----------------------------------------------------------------------
! distribution points are selected; now test for reasonableness
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! has a non-zero distribution area been selected?
!-----------------------------------------------------------------------
if (npts <= 0) then
call document ('init_ms_balance', 'marginal sea ' /&
&/ trim(region_info(n)%name) )
call document ('init_ms_balance', &
'no points selected for distribution area')
call exit_POP (sigAbort, &
'must select at least one set of active points')
endif
!-----------------------------------------------------------------------
! might there be other points within the distribution area
! that were not selected because the search loop exited
! upon reaching the maximum number of points limit?
!-----------------------------------------------------------------------
if (npts == max_points) then
call document ('init_ms_balance', 'marginal sea ' /&
&/ trim(region_info(n)%name) )
call document ('init_ms_balance', &
'warning: an increase of max_points may be necessary')
endif
!-----------------------------------------------------------------------
! are all of the distribution points within one region?
!-----------------------------------------------------------------------
if (my_task == master_task) then
region_number = REGION_MASK_G(ipts(1),jpts(1))
do pt=1,npts
if(region_number /= REGION_MASK_G(ipts(pt),jpts(pt)))then
write(stdout,1000) '(init_ms_balance) ','marginal sea '/&
&/ trim(region_info(n)%name)
write(stdout,1000) '(init_ms_balance)', &
' WARNING: distribution points span two regions'
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
if (pt > 1) then
write(stdout,*)'(init_ms_balance)', ' region1 : ' ,&
REGION_MASK_G(ipts(pt-1),jpts(pt-1))
write(stdout,*) '(init_ms_balance)', ' region2 : ', &
REGION_MASK_G(ipts(pt),jpts(pt))
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
endif
enddo
endif
!-----------------------------------------------------------------------
! define distribution fractions; ensure sum = c1
!-----------------------------------------------------------------------
sum_frac = c0
do pt=1,npts
fracs(pt) = areas(pt)/area
sum_frac = sum_frac + fracs(pt)
enddo
fracs(npts) = fracs(npts) + (c1 - sum_frac)
!-----------------------------------------------------------------------
! define the distribution mask (globally, then locally)
!-----------------------------------------------------------------------
if (my_task == master_task) then
do pt=1,npts
j_loop2: do j=1,ny_global
do i=1,nx_global
if (i == ipts(pt) .and. j == jpts(pt) ) then
MASK_G(i,j) = fracs(pt)
exit j_loop2
endif
enddo
end do j_loop2
enddo
endif
call scatter_global(MASK_FRAC(:,:,:,mask_index),MASK_G,master_task, &
distrb_clinic, &
field_loc_center, field_type_scalar)
region_info(n)%ms_bal%area = area
region_info(n)%ms_bal%mask_index = mask_index
!-----------------------------------------------------------------------
! document the distribution points
!-----------------------------------------------------------------------
if (debug1) then
if (my_task == master_task) then
write(stdout,1002)'(init_ms_balance) ','actual area = ',area
write(stdout,*) trim(region_info(n)%name)
write(stdout,*)'pt ipts jpts TLAT TLON REGION_MASK frac'
do pt=1,npts
write(stdout,*) pt, ipts(pt), jpts(pt), &
TLAT_G(ipts(pt),jpts(pt)), &
TLON_G(ipts(pt),jpts(pt)), &
REGION_MASK_G(ipts(pt),jpts(pt)), &
fracs(pt)
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
enddo
endif
call document('init_ms_balance', 'area = ', area )
call document('init_ms_balance', 'npts = ', npts )
call document('init_ms_balance', 'mask_index = ', mask_index )
endif
!-----------------------------------------------------------------------
! end region_info(n)%marginal_sea if-block
!-----------------------------------------------------------------------
end if
!-----------------------------------------------------------------------
! end marginal-seas do-loop
!-----------------------------------------------------------------------
end do marginal_seas
deallocate (REGION_MASK_G)
deallocate (TLON_G, &
TLAT_G, &
AREAT_G, &
MASK_G )
1000 format(5x,'(',a,')', a )
1002 format(5x,'(',a,') ', a ,1x, 1pe15.5)
!-----------------------------------------------------------------------
!EOC
end subroutine init_ms_balance
!***********************************************************************
!BOP
! !IROUTINE: ms_balancing
! !INTERFACE:
subroutine ms_balancing (STF2, EVAP_F, PREC_F, MELT_F,ROFF_F, IOFF_F, & 2,4
SALT_F, QFLUX, flux_type, ICEOCN_F)
! !DESCRIPTION:
!
! The total excess or deficit of freshwater (kg/s) for each marginal sea
! is transported to or from its associated active-ocean region
!
! The transport term, T, is computed in kg/s of freshwater:
!
!EOP
! T = [ max(0,QFLUX)*c_q +(E+P+M+R+I)*c_f +S*c_s ]*DXT*DYT*c_a
! T == 0 over marginal seas
!
! c_f converts E+P+M+R+I kg/m^2/s freshwater to kg/m^2/s freshwater
! c_f = c1
!
! c_s converts S kg/m^2/s salt to kg/m^2/s freshwater
! c_s = -(1.0e3_r8/ocn_ref_salinity)*rho_fw/rho_sw
! salinity = salt(g)/saltwater(kg)
!
! c_q converts QFLUX W/m^2 heat to kg/m^2/s freshwater
! c_q = -(1.0e4_r8/latent_heat_fusion)*
! (ocn_ref_salinity - sea_ice_salinity)/ocn_ref_salinity
!
! c_a converts T kg/m^2/s * cm^2 to kg/s freshwater
! c_a = 1.0e-4_r8
!
! The total transport is the sum of T kg/s of freshwater. The transport
! is distributed over the designated active-ocean regions
! in freshwater kg/s/cm^2, and must also be converted to the same
! units as STF(:,:,2)
!
! c_t converts flux kg/s/cm^2 to the same units as STF(:,:,2)
! if S(:,:,2) is a salt flux in msu*cm/s
! c_t = -ocn_ref_salinity/rho_fw
! if S(:,:,2) is a freshwater flux in g/cm^2/s
! c_t = 1.0e3_r8
!-----------------------------------------------------------------------
!BOP
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic),intent(in) :: &
EVAP_F ,& ! evaporation flux kg/m^2/s fw
PREC_F ,& ! precipitation flux kg/m^2/s fw
MELT_F ,& ! snow&ice melt flux kg/m^2/s fw
ROFF_F ,& ! river runoff flux kg/m^2/s fw
IOFF_F ,& ! ice runoff flux kg/m^2/s fw
SALT_F ,& ! salt flux kg/m^2/s salt
QFLUX
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), optional, &
intent(in) :: ICEOCN_F
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic),intent(inout):: &
STF2 ! contains STF(:,:,2) from forcing_coupled
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
real (r8) :: &
c_q, &
c_f, &
c_s, &
c_a, &
c_t
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
WORK, &! temporary work array
TRANSPORT ! pointwise transport term
character (*) :: &
flux_type ! flux type of STF(:,:,2) (freshwater or salt )
integer (int_kind) :: &
n, &
iblock
logical (log_kind) :: &
print_transport_daily = .false. ,&
print_transport_monthly = .true. ,&
print_transport_annually = .false. ,&
iceocn_f_present
iceocn_f_present = present(ICEOCN_F)
c_q = -(1.0e4_r8/latent_heat_fusion)* &
(ocn_ref_salinity - sea_ice_salinity)/ocn_ref_salinity
c_f = c1
c_s = -(1.0e3_r8/ocn_ref_salinity)*(rho_fw/rho_sw)
c_a = 1.0e-4_r8
select case (flux_type)
case ('salt')
c_t = -ocn_ref_salinity
case ('freshwater')
c_t = 1.0e3_r8
case default
call exit_POP (sigAbort, &
'ms_balancing -- must set flux_type to either salt or freshwater')
end select
if (eoy) annual_depth = c0
if (eom) monthly_depth = c0
WORK = EVAP_F + PREC_F + MELT_F + ROFF_F + IOFF_F
if ( iceocn_f_present ) WORK = WORK + ICEOCN_F
do n=1,num_regions
TRANSPORT = c0
!-----------------------------------------------------------------------
! in each marginal sea, determine transport term and distribute to
! or from the associated active-ocean region
!-----------------------------------------------------------------------
if (region_info(n)%marginal_sea) then
!-----------------------------------------------------------
! accumulate total transport, in kg/s, for each marginal sea
!
! T = [( max(0,QFLUX)*c_q + (E+P+M+R)*c_f +S*c_s )]*area*c_a
!
! and then set STF2 to zero or the opposite of ice formation
! flux there
!-----------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock = 1,nblocks_clinic
where (REGION_MASK(:,:,iblock) == region_info(n)%number)
TRANSPORT(:,:,iblock) = &
( max(c0,QFLUX(:,:,iblock))*c_q &
+ WORK(:,:,iblock)*c_f &
+ SALT_F(:,:,iblock)*c_s )* &
DXT(:,:,iblock)*DYT(:,:,iblock)*c_a
STF2(:,:,iblock) = -max(c0,QFLUX(:,:,iblock))*c_q*c_t*1.0e-4_r8
end where
enddo ! iblock
!$OMP END PARALLEL DO
region_info(n)%ms_bal%transport = &
global_sum(TRANSPORT,distrb_clinic,field_loc_center)
!-----------------------------------------------------------
! accumulate the average depth of freshwater transported
!-----------------------------------------------------------
annual_depth(n) = annual_depth(n) + &
(region_info(n)%ms_bal%transport/region_info(n)%ms_bal%area) &
*1.0e3_r8*seconds_in_day
monthly_depth(n) = monthly_depth(n) + &
(region_info(n)%ms_bal%transport/region_info(n)%ms_bal%area) &
*1.0e3_r8*seconds_in_day
!----------------------------------------------------------------
! transport excess/deficit to/from associated active-ocean region
!----------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock = 1,nblocks_clinic
STF2(:,:,iblock) = STF2(:,:,iblock) &
+ MASK_FRAC(:,:,iblock,region_info(n)%ms_bal%mask_index) &
* region_info(n)%ms_bal%transport*c_t/ &
(DXT(:,:,iblock)*DYT(:,:,iblock))
enddo ! iblock
!$OMP END PARALLEL DO
endif
enddo ! num_regions
!-----------------------------------------------------------
! record the average daily depth of freshwater transported
! (total transport/total distribution area)
!-----------------------------------------------------------
if (print_transport_daily .and. my_task == master_task) then
write(stdout,*) ' '
call int_to_char (4,iyear ,cyear )
call int_to_char (2,imonth,cmonth)
call int_to_char (2,iday ,cday )
write(stdout,*) cyear/&
&/ '/' /&
&/cmonth/&
&/ '/' /&
&/cday
do n=1,num_regions
if (region_info(n)%marginal_sea) then
write(stdout,1100) &
'average depth of freshwater transported to ', &
trim(region_info(n)%name) /&
&/ ': ', &
(region_info(n)%ms_bal%transport/region_info(n)%ms_bal%area)* &
1.0e3_r8*seconds_in_day, ' (cm/day)'
endif
enddo
write(stdout,*) ' '
endif
!-----------------------------------------------------------
! record the total monthly depth of freshwater transported
!-----------------------------------------------------------
if (print_transport_monthly .and. eom_next &
.and. my_task == master_task) then
write(stdout,*) ' '
do n=1,num_regions
if (region_info(n)%marginal_sea) then
write(stdout,1102) month3_all(imonth) /&
&/ ' ' /&
&/ &
'total monthly depth of freshwater transported to ' , &
trim(region_info(n)%name) /&
&/ ': ',monthly_depth(n),' (cm)'
endif
enddo
write(stdout,*) ' '
endif
!-----------------------------------------------------------
! record the total annual depth of freshwater transported
!-----------------------------------------------------------
if (print_transport_annually .and. eom_next .and. imonth_next == 1 &
.and. my_task == master_task) then
write(stdout,*) ' '
do n=1,num_regions
if (region_info(n)%marginal_sea) then
write(stdout,1101) &
iyear, ' total annual depth of freshwater transported to ' , &
trim(region_info(n)%name) /&
&/ ': ',annual_depth(n),' (cm)'
endif
enddo
write(stdout,*) ' '
endif
1100 format (1x, a45, a20, 1pe25.15, a)
1101 format (1x, i5, a45, a20, 1pe25.15, a)
1102 format (1x, a55, a20, 1pe25.15, a)
!-----------------------------------------------------------------------
!EOC
end subroutine ms_balancing
!***********************************************************************
!BOP
! !IROUTINE: pt_print
! !INTERFACE:
subroutine pt_print (string, i, j) 3
character (*) :: string
integer (int_kind) :: i,j
write(stdout,1000) trim(string), i,j
1000 format (1x,a, 1x, 2i4)
!-----------------------------------------------------------------------
!EOC
end subroutine pt_print
!***********************************************************************
!BOP
! !IROUTINE: print_regions
! !INTERFACE:
subroutine print_regions (REGION_MASK_G) 1
! !DESCRIPTION:
! Print global REGION_MASK
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), dimension(nx_global,ny_global) :: REGION_MASK_G
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
character (4), dimension(25) :: line
integer (int_kind) :: i1, i2, i,j, nn
i1 = 1
i2 = 25
if (my_task == master_task) then
do nn=1,4
write(stdout,2000) 'j ', (i, i=i1,i2)
do j=ny_global,1,-1
do i=i1,i2
if (REGION_MASK_G(i,j) == 0) then
write(line(i-i1+1),'(a4)') '-'
else
write(line(i-i1+1),'(i4)') REGION_MASK_G(i,j)
endif
enddo ! i
write (stdout,2001) j, (line(i-i1+1), i=i1,i2)
enddo ! j
i1 = i2+1; i2 = i1+24
write(stdout,*) ' '
write(stdout,*) ' '
enddo
endif
2001 format(1x, i4,2x, 30a4)
2000 format(2x, a3,2x, 30i4)
!-----------------------------------------------------------------------
!EOC
end subroutine print_regions
end module ms_balance