! subroutines in this file handle the surface mask dependent interpolations
! int22_msk_glb : global area, interpolate from 2-d to 2-d
! int21_msk_glb : global area, interpolate from 2-d to one point
! int22_msk_sub : sub-domain, interpolate from 2-d to 2-d
! int21_msk_sub : sub-domain, interpolate from 2-d to one point
! int2_msk_glb_prep : global area, expand the area (grids) with a specified surface type
! int2_msk_sub_prep : global area, expand the area (grids) with a specified surface type
! dtzm_point : get the vertical mean of the departure from Tf for NSST T-Profile at a point
! dtzm_2d : get the vertical mean of the departure from Tf for NSST T-Profile 2d array
!
subroutine int22_msk_glb(a,isli_a,rlats_a,rlons_a,nlat_a,nlon_a, &
b,isli_b,rlats_b,rlons_b,nlat_b,nlon_b,istyp)
!$$$ subprogram documentation block
! . . .
! subprogram: int22_msk_glb ---(global 2-d array)
! interpolates from a to b with ancillary surface mask (e.g.,
! analysis grid => surface grid) for global arrays
! to gurantee the interpolated value (b) is determined by the
! candidates (a) with the identical surface type from (a)
!
! prgrmmr: li - initial version; org: np2. 02/01/2014
!
! abstract : This routine interpolates a grid to b grid with surface mask accounted
! notes : (1) Here is a 2-d to 2-d interpolation
! (2) The interpolation is performed for specified surface types (istyp)
!
! program history log:
!
! input argument list:
! a - real: 2-d array such as analysis increment at analysis grids
! isli_a - integer: 2-d array: surface mask (0 = water, 1 = land, 2 = sea ice) for a grids
! rlats_a - real: 1-d array: the latitudes of a
! rlons_a - real: 1-d array: the logitudes of a
! nlat_a - integer: number of latitude of a
! nlon_a - integer: number of longitude of a
! isli_b - integer: 2-d array: Analysis surface mask (0 = water, 1 = land, 2 = sea ice) for b grids
! rlats_b - real: 1-d array: the latitudes of b
! rlons_b - real: 1-d array: the logitudes of b
! nlat_b - integer: number of latitude of b
! nlon_b - integer: number of longitude of b
! istyp - integer: target surface type value
!
! output argument list:
! b - real: 2-d array such as analysis increment at surface grids
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block
! USES:
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero,one,two
implicit none
! INPUT:
real (r_kind), dimension(nlat_a,nlon_a), intent(in ) :: a
integer(i_kind), dimension(nlat_a,nlon_a), intent(in ) :: isli_a
real (r_kind), dimension(nlat_a ), intent(in ) :: rlats_a
real (r_kind), dimension(nlon_a ), intent(in ) :: rlons_a
integer(i_kind), dimension(nlat_b,nlon_b), intent(in ) :: isli_b
real (r_kind), dimension(nlat_b ), intent(in ) :: rlats_b
real (r_kind), dimension(nlon_b ), intent(in ) :: rlons_b
integer(i_kind), intent(in ) :: nlat_a,nlon_a,nlat_b,nlon_b,istyp
!OUTPUT:
real (r_kind), dimension(nlat_b,nlon_b), intent( out) :: b
!Declare local variables
integer(i_kind) :: i,j,ix,iy,ii,jj,ixa,iya,sfctyp_b
integer(i_kind) :: nwsum,nfinal
real(r_kind) :: dx0,dx1,dx2,dx3,dy0,dy1,dy2,dy3,dx,dy,dr
real(r_kind) :: dx0s,dx1s,dx2s,dx3s,dy0s,dy1s,dy2s,dy3s
real(r_kind) :: ds00,ds01,ds02,ds03,ds10,ds11,ds12,ds13,ds20,ds21,ds22,ds23,ds30,ds31,ds32,ds33
real(r_kind) :: bavg,bout,dlat,dlon,wsum4,wsum16
real(r_kind), dimension(0:1,0:1) :: w4
real(r_kind), dimension(0:3,0:3) :: w16
dr = 8.0_r_kind ! square of the search radius for 16-point cressman-type analysis
nwsum = 0
nfinal = 0
b=zero
!Loop over all grids of array b to get interpolated value
do j = 1, nlon_b
do i = 1, nlat_b
sfctyp_b = istyp
if ( isli_b(i,j) == sfctyp_b ) then
dlon = rlons_b(j)
call grdcrd1(dlon,rlons_a,nlon_a,1)
iy = int(dlon); dy = dlon-iy; dy1 = one-dy
iy = min(max(0,iy),nlon_a); if(iy == 0) iy = nlon_a
dlat = rlats_b(i)
call grdcrd1(dlat,rlats_a,nlat_a,1)
ix = int(dlat); dx = dlat-ix; dx1 = one-dx
ix = min(max(1,ix),nlat_a)
w4(0,0) = dx1*dy1; w4(0,1) = dx1*dy; w4(1,0) = dx*dy1; w4(1,1) = dx*dy
!
! get the interpolated value with the nearby 4-grids (in a) which has
! the identical surface mask (in b) only
wsum4 = zero
bavg = zero
bout = zero
do jj = 0, 1
iya = iy + jj
if ( iya == nlon_a + 1 ) iya = 1
do ii = 0, 1
ixa = min(nlat_a,ix + ii)
bavg = bavg + w4(ii,jj)*a(ixa,iya)
if ( isli_a(ixa,iya) == sfctyp_b ) then
wsum4 = wsum4 + w4(ii,jj)
bout = bout + w4(ii,jj)*a(ixa,iya)
endif
enddo
enddo
if ( wsum4 > zero ) then
bout = bout/wsum4
else
nwsum = nwsum + 1
! use more candidates from a (extending one more grid futher in both x and y
! direction, which means 16 grids from a will be used) when no same
! surface type candidates can be found in the nearby 4 grids
! to perform a Cressman_type Analysis
ix = ix -1; if(ix == 0) ix = 1
iy = iy -1; if(iy == 0) iy = nlon_a
dx0 = dx + one; dx1 = dx; dx2 = one - dx; dx3 = two - dx
dy0 = dy + one; dy1 = dy; dy2 = one - dy; dy3 = two - dy
dx0s = dx0*dx0; dx1s = dx1*dx1; dx2s = dx2*dx2; dx3s = dx3*dx3
dy0s = dy0*dy0; dy1s = dy1*dy1; dy2s = dy2*dy2; dy3s = dy3*dy3
ds00 = dx0s + dy0s; ds01 = dx0s + dy1s; ds02 = dx0s + dy2s; ds03 = dx0s + dy3s
ds10 = dx1s + dy0s; ds11 = dx1s + dy1s; ds12 = dx1s + dy2s; ds13 = dx1s + dy3s
ds20 = dx2s + dy0s; ds21 = dx2s + dy1s; ds22 = dx2s + dy2s; ds23 = dx2s + dy3s
ds30 = dx3s + dy0s; ds31 = dx3s + dy1s; ds32 = dx3s + dy2s; ds33 = dx3s + dy3s
w16(0,0) = (dr - ds00)/(dr + ds00)
w16(0,1) = (dr - ds01)/(dr + ds01)
w16(0,2) = (dr - ds02)/(dr + ds02)
w16(0,3) = (dr - ds03)/(dr + ds03)
w16(1,0) = (dr - ds10)/(dr + ds10)
w16(1,1) = (dr - ds11)/(dr + ds11)
w16(1,2) = (dr - ds12)/(dr + ds12)
w16(1,3) = (dr - ds13)/(dr + ds13)
w16(2,0) = (dr - ds20)/(dr + ds20)
w16(2,1) = (dr - ds21)/(dr + ds21)
w16(2,2) = (dr - ds22)/(dr + ds22)
w16(2,3) = (dr - ds23)/(dr + ds23)
w16(3,0) = (dr - ds30)/(dr + ds30)
w16(3,1) = (dr - ds31)/(dr + ds31)
w16(3,2) = (dr - ds32)/(dr + ds32)
w16(3,3) = (dr - ds33)/(dr + ds33)
wsum16 = zero
do jj = 0, 3
iya = iy + jj
if ( iya == nlon_a + 1 ) iya = 1
if ( iya == nlon_a + 2 ) iya = 2
if ( iya == nlon_a + 3 ) iya = 3
do ii = 0, 3
ixa = min(nlat_a,ix + ii)
if ( isli_a(ixa,iya) == sfctyp_b ) then
wsum16 = wsum16 + w16(ii,jj)
bout = bout + w16(ii,jj)*a(ixa,iya)
endif
enddo
enddo
if ( wsum16 > zero ) then
bout = bout/wsum16
else
nfinal = nfinal + 1
endif
endif ! if ( wsum4 > zero )
b(i,j)=bout
endif ! if ( isli_b(i,j) == sfctyp_b ) then
enddo ! do i = 1, nlat_b
enddo ! do j = 1, nlon_b
if ( nwsum > 0 ) then
write(*,'(a,I4,I3,I5)') 'int22_msk_glb: Number of grids without specified adjacent surface type, istyp,nwsum ; ',istyp,nwsum
endif
if ( nfinal > 0 ) then
write(*,'(a,I4,I3,I5)') 'int22_msk_glb: Number of grids without interpolted value, istyp,nfinal ; ',istyp,nfinal
endif
end subroutine int22_msk_glb
subroutine int21_msk_sub(a,isli,x,dlat,dlon,istyp,nwsum,nfinal,mype)
!$$$ subprogram documentation block
! . . .
! subprogram: intrp21_msk_sub ---(from 2-d array to obs. location)
!
! prgrmmr: li - initial version; org: np2. 03/01/2014
!
! abstract : This routine interpolates a (2-d array of a subdomain) to a single point with surface mask accounted
! notes : (1) Here is a 2-d to one point interpolation
! (2) The mask is availabe for both 2-d array and the single point
!
! program history log:
!
! input argument list:
! a - real: 2-d array such as analysis increment at analysis grids of a subdomain
! isli - integer: 2-d array: surface mask (0 = water, 1 = land, 2 = sea ice) for grid of a
! dlat - grid relative latitude (obs location)
! dlon - grid relative longitude (obs location)
! istyp - integer: surface type of point x
! mype - mpi task id
! output argument list:
! x - real: a variable (same type of a) at a single point
!$$$ end documentation block
! USES:
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero,one,two
use gridmod, only: istart,jstart,nlon,nlat,lon1,lon2,lat2
implicit none
! INPUT:
real (r_kind), dimension(lat2,lon2), intent(in ) :: a
integer(i_kind), dimension(lat2,lon2), intent(in ) :: isli
real (r_kind), intent(in ) :: dlat,dlon
integer(i_kind), intent(in ) :: istyp,mype
integer(i_kind), intent(inout) :: nwsum,nfinal
!OUTPUT:
real (r_kind), intent( out) :: x
!Declare local variables
integer(i_kind) :: ix1,iy1,ix,iy,ii,jj,ixa,iya,mm1
real(r_kind) :: dx0,dx1,dx2,dx3,dy0,dy1,dy2,dy3,dx,dy,dr
real(r_kind) :: dx0s,dx1s,dx2s,dx3s,dy0s,dy1s,dy2s,dy3s
real(r_kind) :: ds00,ds01,ds02,ds03,ds10,ds11,ds12,ds13,ds20,ds21,ds22,ds23,ds30,ds31,ds32,ds33
real(r_kind) :: bavg,bout,wsum4,wsum16
real(r_kind), dimension(0:1,0:1) :: w4
real(r_kind), dimension(0:3,0:3) :: w16
mm1 = mype + 1
dr = 8.0_r_kind ! square of the search radius for 16-point cressman-type analysis
x=zero
!
! to get interpolated value of x with a (array) and mask info
!
iy1 = int(dlon)
dy = dlon-iy1; dy1 = one-dy
iy =iy1 - jstart(mm1) + 2
if ( iy < 1 ) then
iy1 = iy1 + nlon
iy = iy1 - jstart(mm1) + 2
endif
if ( iy > lon1 + 1 ) then
iy1 = iy1 - nlon
iy = iy1 - jstart(mm1) + 2
endif
ix1 = int(dlat)
ix1 = max(1,min(ix1,nlat))
dx = dlat-ix1; dx1 = one-dx
ix = ix1 - istart(mm1) + 2
w4(0,0) = dx1*dy1; w4(0,1) = dx1*dy; w4(1,0) = dx*dy1; w4(1,1) = dx*dy
!
! get the interpolated value with the nearby 4-grids (in a) which has
! the identical surface mask (istyp for x) only
!
wsum4 = zero
bavg = zero
bout = zero
do jj = 0, 1
iya = min(lon2,iy+jj)
do ii = 0, 1
ixa = min(lat2,ix+ii)
bavg = bavg + w4(ii,jj)*a(ixa,iya)
if ( isli(ixa,iya) == istyp ) then
wsum4 = wsum4 + w4(ii,jj)
bout = bout + w4(ii,jj)*a(ixa,iya)
endif
enddo
enddo
if ( wsum4 > zero ) then
bout = bout/wsum4
else
nwsum = nwsum + 1
! use more candidates from a (extending one more grid futher in both x and y
! direction, which means 16 grids from a will be used) when no same
! surface type candidates can be found in the nearby 4 grids
! to perform a Cressman_type Analysis
ix = ix -1; if(ix == 0) ix = 1
iy = iy -1; if(iy == 0) iy = 1
dx0 = dx + one; dx1 = dx; dx2 = one - dx; dx3 = two - dx
dy0 = dy + one; dy1 = dy; dy2 = one - dy; dy3 = two - dy
dx0s = dx0*dx0; dx1s = dx1*dx1; dx2s = dx2*dx2; dx3s = dx3*dx3
dy0s = dy0*dy0; dy1s = dy1*dy1; dy2s = dy2*dy2; dy3s = dy3*dy3
ds00 = dx0s + dy0s; ds01 = dx0s + dy1s; ds02 = dx0s + dy2s; ds03 = dx0s + dy3s
ds10 = dx1s + dy0s; ds11 = dx1s + dy1s; ds12 = dx1s + dy2s; ds13 = dx1s + dy3s
ds20 = dx2s + dy0s; ds21 = dx2s + dy1s; ds22 = dx2s + dy2s; ds23 = dx2s + dy3s
ds30 = dx3s + dy0s; ds31 = dx3s + dy1s; ds32 = dx3s + dy2s; ds33 = dx3s + dy3s
w16(0,0) = (dr - ds00)/(dr + ds00)
w16(0,1) = (dr - ds01)/(dr + ds01)
w16(0,2) = (dr - ds02)/(dr + ds02)
w16(0,3) = (dr - ds03)/(dr + ds03)
w16(1,0) = (dr - ds10)/(dr + ds10)
w16(1,1) = (dr - ds11)/(dr + ds11)
w16(1,2) = (dr - ds12)/(dr + ds12)
w16(1,3) = (dr - ds13)/(dr + ds13)
w16(2,0) = (dr - ds20)/(dr + ds20)
w16(2,1) = (dr - ds21)/(dr + ds21)
w16(2,2) = (dr - ds22)/(dr + ds22)
w16(2,3) = (dr - ds23)/(dr + ds23)
w16(3,0) = (dr - ds30)/(dr + ds30)
w16(3,1) = (dr - ds31)/(dr + ds31)
w16(3,2) = (dr - ds32)/(dr + ds32)
w16(3,3) = (dr - ds33)/(dr + ds33)
wsum16 = zero
do jj = 0, 3
iya = min(lon2,iy+jj)
do ii = 0, 3
ixa = min(lat2,ix + ii)
if ( isli(ixa,iya) == istyp ) then
wsum16 = wsum16 + w16(ii,jj)
bout = bout + w16(ii,jj)*a(ixa,iya)
endif
enddo
enddo
if ( wsum16 > zero ) then
bout = bout/wsum16
else
nfinal = nfinal + 1
endif
endif ! if ( wsum4 > zero )
x=bout
! if ( nwsum > 0 ) then
! write(*,'(a,I4,I3,I5)') 'int21_msk_sub: Number of grids without specified adjacent surface type, mype,istyp,nwsum ; ', mype,istyp,nwsum
! endif
! if ( nfinal > 0 ) then
! write(*,'(a,I4,I3,I5)') 'int21_msk_sub: Number of grids without interpolted value, mype,istyp,nfinal ; ', mype,istyp,nfinal
! endif
end subroutine int21_msk_sub
subroutine int2_msk_glb_prep(a,isli_a,b,isli_b,nx,ny,sfctyp_b,nprep)
!$$$ subprogram documentation block
! . . .
! subroutine: int2_msk_glb_prep --- (global 2-d array)
! for a specified surface type (sfctyp_b), expanding the area (grids) with this
! surface type using the surounded 8 grids with the same surface type
!
! prgrmmr: li - initial version; org: np2
!
! input argument list:
! a - real: 2-d array
! isli_a - integer: 2-d array: surface mask (0 = water, 1 = land, 2 = seaice) for a grids
! nx - integer: number of grids in x-direction
! ny - integer: number of grids in y-direction
! sfctyp_b - integer: the targeted surface type (0=water, 1=land, 2-sea ice)
! nprep - integer: number of times to do the extension
!
! output argument list:
! b - real: 2-d array such as analysis increment at surface grids
! isli_b - integer: 2-d array such as analysis surface mask
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
! USES:
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero,one,two
implicit none
! INPUT:
real (r_kind), dimension(nx,ny), intent(in ) :: a
integer(i_kind), dimension(nx,ny), intent(in ) :: isli_a
integer(i_kind), intent(in ) :: nx,ny,sfctyp_b,nprep
!OUTPUT:
real (r_kind), dimension(nx,ny), intent( out) :: b
integer(i_kind), dimension(nx,ny), intent( out) :: isli_b
!Declare local variables
integer(i_kind) :: i,j,k,ii,jj,ix,iy,n
real(r_kind) :: bout
real(r_kind), dimension(nx,ny) :: wb
integer(r_kind), dimension(nx,ny) :: wmsk
!
! Initialize b and isli_b
!
b = a; isli_b = isli_a
do k = 1, nprep
!
! Initialize/update work array for the value and mask
!
wb = b; wmsk = isli_b
!
! Loop over all grids of array b to update the grids nearby the sfctyp_b surface type
!
do j = 1, ny
do i = 1, nx
if ( wmsk(i,j) /= sfctyp_b ) then
n = 0
bout = zero
do jj = j - 1, j + 1
do ii = i - 1, i + 1
iy = jj; if (iy == ny + 1) iy = 1; if (iy == 0) iy = ny
ix = min(max(ii,1),nx)
if ( wmsk(ix,iy) == sfctyp_b ) then
n = n + 1
bout = bout + wb(ix,iy)
endif
enddo
enddo
if ( n > 0 ) then
b(i,j) = bout/real(n)
isli_b(i,j) = sfctyp_b
endif
endif
enddo
enddo
enddo ! do k = 1, nprep
end subroutine int2_msk_glb_prep
subroutine int2_msk_sub_prep(a,isli_a,b,isli_b,nx,ny,sfctyp_b,nprep)
!$$$ subprogram documentation block
! . . .
! subroutine: int2_msk_sub_prep --- (subdomain 2-d array)
! for a specified surface type (sfctyp_b), expanding the area (grids) with this
! surface type using the surounded 8 grids with the same surface type
!
! prgrmmr: li - initial version; org: np2
!
! input argument list:
! a - real: 2-d array
! isli_a - integer: 2-d array: surface mask (0 = water, 1 = land, 2 = seaice) for a grids
! nx - integer: number of grids in x-direction
! ny - integer: number of grids in y-direction
! sfctyp_b - integer: the targeted surface type (0=water, 1=land, 2-sea ice)
! nprep - integer: number of times to do the extension
!
! output argument list:
! b - real: 2-d array such as analysis increment at surface grids
! isli_b - integer: 2-d array such as analysis surface mask
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
! USES:
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero,one,two
implicit none
! INPUT:
real (r_kind), dimension(nx,ny), intent(in ) :: a
integer(i_kind), dimension(nx,ny), intent(in ) :: isli_a
integer(i_kind), intent(in ) :: nx,ny,sfctyp_b,nprep
!OUTPUT:
real (r_kind), dimension(nx,ny), intent( out) :: b
integer(i_kind), dimension(nx,ny), intent( out) :: isli_b
!Declare local variables
integer(i_kind) :: i,j,k,ii,jj,ix,iy,n
real(r_kind) :: bout
real(r_kind), dimension(nx,ny) :: wb
integer(r_kind), dimension(nx,ny) :: wmsk
!
! Initialize b and isli_b
!
b = a; isli_b = isli_a
do k = 1, nprep
!
! Initialize/update work array for the value and mask
!
wb = b; wmsk = isli_b
!
! Loop over all grids of array b to update the grids nearby the sfctyp_b surface type
!
do j = 1, ny
do i = 1, nx
if ( wmsk(i,j) /= sfctyp_b ) then
n = 0
bout = zero
do jj = j - 1, j + 1
do ii = i - 1, i + 1
iy = min(max(jj,1),ny)
ix = min(max(ii,1),nx)
if ( wmsk(ix,iy) == sfctyp_b ) then
n = n + 1
bout = bout + wb(ix,iy)
endif
enddo
enddo
if ( n > 0 ) then
b(i,j) = bout/real(n)
isli_b(i,j) = sfctyp_b
endif
endif
enddo
enddo
enddo ! do k = 1, nprep
end subroutine int2_msk_sub_prep
!*******************************************************************************************
subroutine dtzm_point(xt,xz,dt_cool,zc,z1,z2,dtzm)
! ===================================================================== !
! !
! description: get dtzm = mean of dT(z) (z1 - z2) with NSST dT(z) !
! dT(z) = (1-z/xz)*dt_warm - (1-z/zc)*dt_cool !
! !
! usage: !
! !
! call dtzm_point !
! !
! inputs: !
! (xt,xz,dt_cool,zc,z1,z2, !
! outputs: !
! dtzm) !
! !
! program history log: !
! !
! 2015 -- xu li createad original code !
! inputs: !
! xt - real, heat content in dtl 1 !
! xz - real, dtl thickness 1 !
! dt_cool - real, sub-layer cooling amount 1 !
! zc - sub-layer cooling thickness 1 !
! z1 - lower bound of depth of sea temperature 1 !
! z2 - upper bound of depth of sea temperature 1 !
! outputs: !
! dtzm - mean of dT(z) (z1 to z2) 1 !
!
use kinds, only: r_single, i_kind
use constants, only: zero,one,half
implicit none
real (kind=r_single), intent(in) :: xt,xz,dt_cool,zc,z1,z2
real (kind=r_single), intent(out) :: dtzm
! Local variables
real (kind=r_single) :: dt_warm,dtw,dtc
!
! get the mean warming in the range of z=z1 to z=z2
!
dtw = zero
if ( xt > zero ) then
dt_warm = (xt+xt)/xz ! Tw(0)
if ( z1 < z2) then
if ( z2 < xz ) then
dtw = dt_warm*(one-(z1+z2)/(xz+xz))
elseif ( z1 < xz .and. z2 >= xz ) then
dtw = half*(one-z1/xz)*dt_warm*(xz-z1)/(z2-z1)
endif
elseif ( z1 == z2 ) then
if ( z1 < xz ) then
dtw = dt_warm*(one-z1/xz)
endif
endif
endif
!
! get the mean cooling in the range of z=z1 to z=z2
!
dtc = zero
if ( zc > zero ) then
if ( z1 < z2) then
if ( z2 < zc ) then
dtc = dt_cool*(one-(z1+z2)/(zc+zc))
elseif ( z1 < zc .and. z2 >= zc ) then
dtc = half*(one-z1/zc)*dt_cool*(zc-z1)/(z2-z1)
endif
elseif ( z1 == z2 ) then
if ( z1 < zc ) then
dtc = dt_cool*(one-z1/zc)
endif
endif
endif
!
! get the mean T departure from Tf in the range of z=z1 to z=z2
!
dtzm = dtw - dtc
end subroutine dtzm_point
!*******************************************************************************************
!*******************************************************************************************
subroutine dtzm_2d(xt,xz,dt_cool,zc,slmsk,z1,z2,nx,ny,dtzm)
! ===================================================================== !
! !
! description: get dtzm = mean of dT(z) (z1 - z2) with NSST dT(z) !
! dT(z) = (1-z/xz)*dt_warm - (1-z/zc)*dt_cool !
! !
! usage: !
! !
! call dtzm_2d !
! !
! inputs: !
! (xt,xz,dt_cool,zc,z1,z2, !
! outputs: !
! dtzm) !
! !
! program history log: !
! !
! 2015 -- xu li createad original code !
! inputs: !
! xt - real, heat content in dtl !
! xz - real, dtl thickness !
! dt_cool - real, sub-layer cooling amount !
! zc - sub-layer cooling thickness !
! nx - integer, dimension in x-direction (zonal) !
! ny - integer, dimension in y-direction (meridional) !
! z1 - lower bound of depth of sea temperature !
! z2 - upper bound of depth of sea temperature !
! outputs: !
! dtzm - mean of dT(z) (z1 to z2) !
!
use kinds, only: r_single, i_kind
use constants, only: zero,one,half
implicit none
real(kind=r_single), dimension(nx,ny), intent(in) :: xt,xz,dt_cool,zc,slmsk
real(kind=r_single), intent(in) :: z1,z2
integer(kind=i_kind), intent(in) :: nx,ny
real(kind=r_single), dimension(nx,ny), intent(out) :: dtzm
! Local variables
integer(kind=i_kind) :: i,j
real(kind=r_single), dimension(nx,ny) :: dtw,dtc
real(kind=r_single) :: dt_warm
!
! initialize dtw & dtc as zeros
!
dtw(:,:) = zero
dtc(:,:) = zero
do j = 1, ny
do i= 1, nx
if ( slmsk(i,j) == zero ) then
!
! get the mean warming in the range of z=z1 to z=z2
!
if ( xt(i,j) > zero ) then
dt_warm = (xt(i,j)+xt(i,j))/xz(i,j) ! Tw(0)
if ( z1 < z2) then
if ( z2 < xz(i,j) ) then
dtw(i,j) = dt_warm*(one-(z1+z2)/(xz(i,j)+xz(i,j)))
elseif ( z1 < xz(i,j) .and. z2 >= xz(i,j) ) then
dtw(i,j) = half*(one-z1/xz(i,j))*dt_warm*(xz(i,j)-z1)/(z2-z1)
endif
elseif ( z1 == z2 ) then
if ( z1 < xz(i,j) ) then
dtw(i,j) = dt_warm*(one-z1/xz(i,j))
endif
endif
endif
!
! get the mean cooling in the range of z=0 to z=zsea
!
if ( zc(i,j) > zero ) then
if ( z1 < z2) then
if ( z2 < zc(i,j) ) then
dtc(i,j) = dt_cool(i,j)*(one-(z1+z2)/(zc(i,j)+zc(i,j)))
elseif ( z1 < zc(i,j) .and. z2 >= zc(i,j) ) then
dtc(i,j) = half*(one-z1/zc(i,j))*dt_cool(i,j)*(zc(i,j)-z1)/(z2-z1)
endif
elseif ( z1 == z2 ) then
if ( z1 < zc(i,j) ) then
dtc(i,j) = dt_cool(i,j)*(one-z1/zc(i,j))
endif
endif
endif
endif ! if ( slmsk(i,j) == zero ) then
enddo
enddo
!
! get the mean T departure from Tf in the range of z=z1 to z=z2
!
where ( slmsk(:,:) == zero )
dtzm(:,:) = dtw(:,:) - dtc(:,:)
end where
end subroutine dtzm_2d