module gsd_update_mod
!$$$ module documentation block
! . . . .
! module: gsd_update_mod module for updating surface, soil, moisture from GSD for RR
! prgmmr: Hu org: gsd date: 2012-01-12
!
! abstract: module for updating surface, soil, moisture from GSD for RR
!
! program history log:
! 2012-01-12 Hu
! 2015-01-12 Hu fix the bug in coast proximity calculation in subdomain
! 2015-01-14 Hu do T soil nudging over snow
! 2015-01-15 Hu move the land/sea mask check to fine grid update step
!
! subroutines included:
! sub gsd_update_soil_tq - change surface and soil based on analysis increment
! sub gsd_limit_ocean_q - limits to analysis increments over oceans
! sub gsd_update_th2 - adjust 2-m t based on analysis increment
! sub gsd_update_q2 - adjust 2-m q based on analysis increment
!
! Variable Definitions:
use gsi_metguess_mod, only: GSI_MetGuess_Bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use mpeu_util, only: die
implicit none
! set default to private
private
! set subroutines to public
public :: gsd_update_soil_tq
public :: gsd_limit_ocean_q
public :: gsd_update_th2
public :: gsd_update_q2
public :: gsd_gen_coast_prox
! set passed variables to public
contains
subroutine gsd_update_soil_tq(tinc,is_t,qinc,is_q,it)
!$$$ subprogram documentation block
! . . . .
! subprogram: update_surface change surface and soil based on analysis increment
! prgmmr: Hu org: GSD date: 2011-08-18
!
! abstract: This routine does the following things:
! 1) add lowest level t increment to T2
!
!
! program history log:
! 1990-10-06 parrish - original code
! 2013-10-19 todling - metguess now holds background
!
! input argument list:
! tinc : first level temperature analysis increment
! qinc : first level moisture analysis increment
!
! output argument list:
!
! comments:
!
! attributes:
!$$$
use kinds, only: r_kind,i_kind
use jfunc, only: tsensible,qoption
use derivsmod, only: qsatg
use constants, only: zero,one,fv,one_tenth,deg2rad,pi
use constants, only: partialSnowThreshold,t0c
use gridmod, only: lat2,lon2,nsig,nsig_soil
use gridmod, only: regional_time
use guess_grids, only: ges_tsen,sno,coast_prox
use wrf_mass_guess_mod, only: ges_xlon,ges_xlat
use guess_grids, only: ges_prsl
! use guess_grids, only: nfldsig
use rapidrefresh_cldsurf_mod, only: l_gsd_soiltq_nudge
implicit none
! Declare passed variables
integer(i_kind) is_t,is_q
real(r_kind),dimension(lat2,lon2), intent(in) :: tinc
real(r_kind),dimension(lat2,lon2), intent(in) :: qinc
integer,intent(in) :: it ! guess time level
! Declare local variables
real(r_kind),dimension(lat2,lon2) :: csza
integer(i_kind) :: gmt,nday,iyear,imonth,iday
real(r_kind) :: declin
real(r_kind) :: hrang,xxlat
real(r_kind) :: sumqc
logical ice
integer(i_kind) :: iderivative
real(r_kind),allocatable,dimension(:,:,:):: rhgues
integer(i_kind) i,j,k,ier,istatus
real(r_kind) :: ainc, tinct
real(r_kind) :: coast_fac,temp,temp_fac,dts_min,tincf
real(r_kind) :: snowthreshold
!
real(r_kind), pointer :: ges_qc(:,:,:) ! cloud water
real(r_kind), pointer :: ges_qi(:,:,:) ! could ice
real(r_kind),dimension(:,: ),pointer:: ges_tsk =>NULL()
real(r_kind),dimension(:,: ),pointer:: ges_soilt1=>NULL()
real(r_kind),dimension(:,:,:),pointer:: ges_tslb =>NULL()
real(r_kind),dimension(:,:,:),pointer:: ges_smois =>NULL()
real(r_kind),dimension(:,:,:),pointer:: ges_q =>NULL()
integer(i_kind) :: itsig
!*******************************************************************************
!
snowthreshold=1.0e-10_r_kind
itsig=it
ier=0
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'ql',ges_qc,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'qi',ges_qi,istatus);ier=ier+istatus
if(ier/=0) then
write(6,*) 'No cloud water and ice for soil nudging!'
return ! no guess, nothing to do
endif
!
! calculation solar declination
!
iyear=regional_time(1)
imonth=regional_time(2)
iday=regional_time(3)
call w3fs13(iyear,imonth,iday,nday)
declin=deg2rad*23.45_r_kind*sin(2.0_r_kind*pi*(284+nday)/365.0_r_kind)
!
! csza = fraction of solar constant (cos of zenith angle)
gmt = regional_time(4) ! UTC
do j=2,lon2
do i=1,lat2
hrang=15._r_kind*gmt*deg2rad + ges_xlon(i,j,1) - pi
xxlat=ges_xlat(i,j,1)
csza(i,j)=sin(xxlat)*sin(declin) &
+cos(xxlat)*cos(declin)*cos(hrang)
end do
end do
if( l_gsd_soiltq_nudge .and. is_t > 0) then
! --------------------------------------------
! --- Increment top level of soil temp and snow temp
! ONLY AT LAND POINTS according to
! sfc temp increment.
! --------------------------------------------
! -- modifications to reintroduce soil temperature nudging
! -- Stan and Tanya - 15 July 2004
! -- allow cooling of soil only (not warming)
! - allow only up to 1.0 K (half of negative ainc for temp)
! - don't allow if snow water is > 6 mm
! do it=1,nfldsig
ier=0
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'tskn' ,ges_tsk ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'tsoil',ges_soilt1,istatus)
ier=ier+istatus
call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'smoist',ges_smois,istatus)
ier=ier+istatus
call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'tslb' ,ges_tslb ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'q' ,ges_q ,istatus)
ier=ier+istatus
if(ier/=0) then
! code doesn't have to die here ... needs attention for generalization
call die('gsd_update_soil_tq',': cannot find q in guess, ier=',ier)
endif
do j=1,lon2
do i=1,lat2
if(tsensible) then
ainc=tinc(i,j)
else
ainc=tinc(i,j)/(one+fv*ges_q(i,j,1))
endif
coast_fac = max(0._r_kind,(coast_prox(i,j)-0.5_r_kind))/0.5_r_kind
temp = ges_tsen(i,j,1,it)
! *** Increase soil adjustment by factor of 2.0 if temps are
! 25 deg C, and 2.5 if temp is 32.5 C .
! -- Stan B. (John, Tanya) - 31 July 2006
temp_fac = 1._r_kind+min (1.5_r_kind,max(0._r_kind,(temp-283._r_kind)/15._r_kind))
dts_min = -2._r_kind
! -- Allow soil temp cooling to be up to 2.5 * 0.6 = 1.5 X 2 C ( = 3K)
dts_min = dts_min*temp_fac*0.6_r_kind
! mhu, Jan 15,2015: move the land/sea masck check to fine grid update step
tincf = ainc*temp_fac*coast_fac
! mhu and Tanya: Jan 14, 2015: do T soil nudging over snow
if(nsig_soil == 9) then
! - top level soil temp
ges_tslb(i,j,1) = ges_tslb(i,j,1) + &
min(1._r_kind,max(dts_min,tincf*0.6_r_kind))
! - 0-1 cm level - soil temp
ges_tslb(i,j,2) = ges_tslb(i,j,2) + &
min(1._r_kind,max(dts_min,tincf*0.55_r_kind))
! - 1-4 cm level - soil temp
ges_tslb(i,j,3) = ges_tslb(i,j,3) + &
min(1._r_kind,max(dts_min,tincf*0.4_r_kind))
! - 4-10 cm level - soil temp
ges_tslb(i,j,4) = ges_tslb(i,j,4) + &
min(1._r_kind,max(dts_min,tincf*0.3_r_kind))
! - 10-30 cm level - soil temp
ges_tslb(i,j,5) = ges_tslb(i,j,5) + &
min(1._r_kind,max(dts_min,tincf*0.2_r_kind))
else
! - top level soil temp
ges_tslb(i,j,1) = ges_tslb(i,j,1) + &
min(1._r_kind,max(dts_min,tincf*0.6_r_kind))
! - 0-5 cm level - soil temp
ges_tslb(i,j,2) = ges_tslb(i,j,2) + &
min(1._r_kind,max(dts_min,tincf*0.4_r_kind))
! - 5-20 cm level - soil temp
ges_tslb(i,j,3) = ges_tslb(i,j,3) + &
min(1._r_kind,max(dts_min,tincf*0.2_r_kind))
endif
if (sno(i,j,it) < partialSnowThreshold) THEN
! partialSnowThreshold (32 mm) is the threshold for partial snow.
! When grid cell is partially covered with snow or snow-free - always update TSK and SOILT1
ges_tsk(i,j) = ges_tsk(i,j) + min(1._r_kind,max(dts_min,tincf*0.6_r_kind))
ges_soilt1(i,j) = ges_soilt1(i,j) + min(1._r_kind,max(dts_min,tincf*0.6_r_kind))
else
! grid cell is fully covered with snow
if(tincf < zero) then
! always adjust TSK and SOILT1 when tincf < 0 - cooling
ges_tsk(i,j) = ges_tsk(i,j) + min(1._r_kind,max(-2._r_kind,tincf*0.6_r_kind))
ges_soilt1(i,j) = ges_soilt1(i,j) + min(1._r_kind,max(-2._r_kind,tincf*0.6_r_kind))
else
! if ticnf > 0 - warming, then adjust snow TSK and SOILT1 only if TSK < t0c (273 K).
! If TSK > t0c(273 K) most likely due to melting process, then leave TSK and SOILT1 unchanged.
if(ges_tsk(i,j) < t0c ) then
ges_tsk(i,j) = min(t0c,ges_tsk(i,j) + min(1._r_kind,max(-2._r_kind,tincf*0.6_r_kind)))
ges_soilt1(i,j) = min(t0c,ges_soilt1(i,j) + min(1._r_kind,max(-2._r_kind,tincf*0.6_r_kind)))
endif ! tsk < 273 K
endif ! tincf < 0.
endif ! sno(i,j,it) < 32
end do
end do
! end do ! it
!---------------------------------------------------------
! Nudge soil moisture
! Tanya S. and Stan B. - 21 July 2004 - first version
!---------------------------------------------------------
! Compute saturation specific humidity.
iderivative = 0
if(qoption == 1 )then
iderivative = 1
else
iderivative = 2
end if
ice=.true.
call genqsat(qsatg,ges_tsen(1,1,1,it),ges_prsl(1,1,1,it), &
lat2,lon2,nsig,ice,iderivative)
allocate(rhgues(lat2,lon2,nsig))
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'q',ges_q,istatus)
if(istatus/=0) then
! code doesn't have to die here ... needs attention for generalization
call die('gsd_update_soil_tq',': cannot find q in guess')
endif
do k=1,nsig
do j=1,lon2
do i=1,lat2
rhgues(i,j,k)=ges_q(i,j,k)/qsatg(i,j,k)
end do
end do
end do
if( is_q > 0) then
! do it=1,nfldsig
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'q',ges_q,istatus)
ier=istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'smoist',ges_smois,istatus)
ier=ier+istatus
if(ier/=0) then
! code doesn't have to die here ... needs attention for generalization
call die('gsd_update_soil_tq',': cannot find q in guess')
endif
do j=1,lon2
do i=1,lat2
if(tsensible) then
tinct=tinc(i,j)
else
tinct=tinc(i,j)/(one+fv*ges_q(i,j,1))
endif
ainc=qinc(i,j)/qsatg(i,j,1) ! analysis increment in RH
! -- use overall limits based on k level
ainc = max(-0.3_r_kind,min(0.3_r_kind,ainc))
! -- When background is already dry and pRH increment
! is negative (toward drier still), limit ainc further.
if (rhgues(i,j,1) < 0.2_r_kind .and. ainc < 0.0_r_kind ) then
ainc=ainc*rhgues(i,j,1)/0.2_r_kind
end if
if (rhgues(i,j,1) < 0.4_r_kind .and. ainc < 0.0_r_kind ) then
ainc=ainc*rhgues(i,j,1)/0.4_r_kind
end if
ainc = max(-0.15_r_kind,min(0.15_r_kind,ainc))
! - Only do nudging over land, if daytime (defined as
! cos of sun zenith angle > 0.1), and if
! sfc temp increment is negative (meaning that
! background sfc temp was too warm)
! -- some adjustments below to soil moisture adjustment,
! which seems to have resulted in too much moistening
! overall. Stan B. - 24 Oct 04 - 04z
! mhu, Jan 15,2015: move the land/sea masck check to fine grid update step
! if (isli(i,j,it) == 1 .and. csza(i,j) > 0.3_r_kind) then
if (csza(i,j) > 0.3_r_kind) then
sumqc=0
do k=1,nsig
sumqc=max(sumqc,max(ges_qc(i,j,k),ges_qi(i,j,k)))
enddo
sumqc=0 ! trun off cloud
if( sumqc < 1.0e-6_r_kind) then
if( sno(i,j,it) < snowthreshold ) then ! don't do the
! moisture adjustment if there is snow
if (tinct < -0.15_r_kind) then
! - top level soil moisture
! -- mod - 3/15/13
! increase moistening from factor of 0.2 to 0.3
ges_smois(i,j,1) = min (max(ges_smois(i,j,1),ges_smois(i,j,2)), &
ges_smois(i,j,1) + min(0.03_r_kind,max(0._r_kind,(ainc*0.2_r_kind))))
ges_smois(i,j,2) = min (max(ges_smois(i,j,2),ges_smois(i,j,3)), &
ges_smois(i,j,2) + min(0.03_r_kind,max(0._r_kind,(ainc*0.2_r_kind))))
if(nsig_soil == 9) then
ges_smois(i,j,3) = min (max(ges_smois(i,j,3),ges_smois(i,j,4)), &
ges_smois(i,j,3) + min(0.03_r_kind,max(0._r_kind,(ainc*0.2_r_kind))))
ges_smois(i,j,4) = min (max(ges_smois(i,j,4),ges_smois(i,j,5)), &
ges_smois(i,j,4) + min(0.03_r_kind,max(0._r_kind,(ainc*0.1_r_kind))))
endif
! -- above logic
! 7/26/04 -
! previously - min was sm1_p (level 2)
! now - min is (max of level 1 and level 2)
! Implication - If level 1 was already more moist than
! level 2, don't force level 1 SM back down to level 2.
! -- mod - 5/1/05
! Decrease moistening from factor of 0.2 to 0.1
! -- mod - 3/15/13
! increase moistening from factor of 0.1 to 0.3
endif
if (tinct > 0.15_r_kind) then
! - top level soil moisture
! -- addition 5/1/05
! Now also dry soil if tinc is positive (warming)
! and the RH_inc is negative.
ges_smois(i,j,1) = max(0.0_r_kind,ges_smois(i,j,1) + &
max(-0.03_r_kind,min(0._r_kind,(ainc*0.2_r_kind))))
ges_smois(i,j,2) = max(0.0_r_kind,ges_smois(i,j,2) + &
max(-0.03_r_kind,min(0._r_kind,(ainc*0.2_r_kind))))
if(nsig_soil == 9) then
ges_smois(i,j,3) = max(0.0_r_kind,ges_smois(i,j,3) + &
max(-0.03_r_kind,min(0._r_kind,(ainc*0.2_r_kind))))
ges_smois(i,j,4) = max(0.0_r_kind,ges_smois(i,j,4) + &
max(-0.03_r_kind,min(0._r_kind,(ainc*0.1_r_kind))))
endif
END IF
endif ! sno(i,j,it) < snowthreshold
endif ! sumqc < 1.0e-6_r_kind
endif
end do
end do
! end do ! it
endif ! is_q > 0
deallocate(rhgues)
endif
return
end subroutine gsd_update_soil_tq
subroutine gsd_limit_ocean_q(qinc,it)
!$$$ subprogram documentation block
! . . . .
! subprogram: gsd_limit_ocean_q Rlimits to analysis increments over oceans
! prgmmr: Hu org: GSD date: 2011-08-31
!
! abstract: This routine does the following things:
!
!
! program history log:
! 2011-08-31 Hu - original code
!
! input argument list:
! qinc : moisture analysis increment
!
! output argument list:
!
! comments:
!
! attributes:
!$$$
use kinds, only: r_kind,i_kind
use jfunc, only: qoption
use derivsmod, only: qsatg
use constants, only: zero,one,one_tenth
use gridmod, only: lat2,lon2,nsig
use guess_grids, only: ges_tsen,ges_prsl
use guess_grids, only: isli
! use guess_grids, only: nfldsig
implicit none
! Declare passed variables
integer(i_kind) istatus
real(r_kind),dimension(lat2,lon2,nsig), intent(inout) :: qinc
integer,intent(in) :: it ! guess time level
! Declare local variables
logical ice
integer(i_kind) :: i,j,k,iderivative
real(r_kind),allocatable,dimension(:,:,:):: rhgues
real(r_kind) :: qinc_rh
real(r_kind),dimension(:,:,:),pointer:: ges_q=>NULL()
! Compute saturation specific humidity.
iderivative = 0
if(qoption == 1 )then
iderivative = 1
else
iderivative = 2
end if
ice=.true.
call genqsat(qsatg,ges_tsen(1,1,1,it),ges_prsl(1,1,1,it),lat2,lon2,nsig,ice,iderivative)
allocate(rhgues(lat2,lon2,nsig))
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(it),'q',ges_q,istatus)
if(istatus/=0) then
! code doesn't have to die here ... needs attention for generalization
call die('gsd_update_soil_tq',': cannot find q in guess')
endif
do k=1,nsig
do j=1,lon2
do i=1,lat2
rhgues(i,j,k)=ges_q(i,j,k)/qsatg(i,j,k)
end do
end do
end do
! do it=1,nfldsig
do k=1,nsig
do j=1,lon2
do i=1,lat2
if (isli(i,j,1) == 0) then
qinc_rh=qinc(i,j,k)/qsatg(i,j,k)
! -- limit pseudo-RH change over water to +/- 0.1
if( k <= 10) then
qinc_rh=max(-0.1_r_kind,min(0.1_r_kind,qinc_rh))
elseif( k >= 11.and.k <=18 ) then
qinc_rh=max(-0.2_r_kind,min(0.2_r_kind,qinc_rh))
endif
! -- Limit further drying out over water and near surface.
if(rhgues(i,j,k) < 0.6_r_kind .and. k <=4 .and. qinc_rh < zero) then
qinc_rh=qinc_rh*rhgues(i,j,k)/1.0_r_kind
endif
qinc(i,j,k)=qinc_rh*qsatg(i,j,k)
else
qinc(i,j,k)=qinc(i,j,k)
end if ! isli(i,j,1)
end do
end do
end do
! end do
deallocate(rhgues)
end subroutine gsd_limit_ocean_q
subroutine gsd_update_th2(tinc,it)
!$$$ subprogram documentation block
! . . . .
! subprogram: gsd_update_th2 adjust 2-m t based on analysis increment
! prgmmr: Hu org: GSD date: 2011-10-04
!
! abstract: This routine does the following things:
! 1) add lowest level t increment to T2
!
!
! program history log:
! 2011-10-04 parrish - original code
! 2013-10-19 todling - get guess fileds from bundle
! 2017-03-23 Hu - add code to use hybrid vertical coodinate in WRF
! MASS core
!
! input argument list:
! tinc : first level temperature analysis increment
!
! output argument list:
!
! comments:
!
! attributes:
!$$$
use kinds, only: r_kind,i_kind
use jfunc, only: tsensible
use constants, only: zero,one,fv,rd_over_cp_mass,one_tenth
use gridmod, only: lat2,lon2,aeta1_ll,pt_ll,aeta2_ll
! use guess_grids, only: nfldsig
implicit none
! Declare passed variables
real(r_kind),dimension(lat2,lon2), intent(in) :: tinc
integer,intent(in) :: it ! guess time level
real(r_kind),parameter:: r10=10.0_r_kind
real(r_kind),parameter:: r100=100.0_r_kind
integer(i_kind) i,j,ier,ihaveq
real(r_kind) :: dth2, work_prsl,work_prslk
real(r_kind),dimension(:,: ),pointer:: ges_ps =>NULL()
real(r_kind),dimension(:,: ),pointer:: ges_th2=>NULL()
real(r_kind),dimension(:,:,:),pointer:: ges_q =>NULL()
!*******************************************************************************
!
! 2-m temperature
! do it=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(it),'th2m',ges_th2,ier)
if(ier/=0) return
call gsi_bundlegetpointer(gsi_metguess_bundle(it),'ps',ges_ps,ier)
if(ier/=0) return
! NOTE: for some odd reason the orig. code before bundle change was getting q
! from slot it=1 - to preserve zero diff I left as such - RTodling
call gsi_bundlegetpointer(gsi_metguess_bundle(it),'q' ,ges_q ,ihaveq)
do j=1,lon2
do i=1,lat2
if(tsensible) then
dth2=tinc(i,j)
else
if(ihaveq/=0) cycle
dth2=tinc(i,j)/(one+fv*ges_q(i,j,1))
endif
! Convert sensible temperature to potential temperature
work_prsl = one_tenth*(aeta1_ll(1)*(r10*ges_ps(i,j)-pt_ll)+ &
aeta2_ll(1) + pt_ll)
work_prslk = (work_prsl/r100)**rd_over_cp_mass
ges_th2(i,j) = ges_th2(i,j) + dth2/work_prslk
end do
end do
! end do
return
end subroutine gsd_update_th2
subroutine gsd_update_q2(qinc,it)
!$$$ subprogram documentation block
! . . . .
! subprogram: gsd_update_q2 adjust 2-m q based on analysis increment
! prgmmr: Hu org: GSD date: 2011-10-04
!
! abstract: This routine does the following things:
! 1) add lowest level q increment to Q2
!
!
! program history log:
! 2014-01-22 Hu - original code
! 2014-04-04 Todling - ges_q2 now in MetGuess
!
! input argument list:
! qinc : first level moisture analysis increment
!
! output argument list:
!
! comments:
!
! attributes:
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2
! use guess_grids, only: nfldsig
implicit none
! Declare passed variables
real(r_kind),dimension(lat2,lon2), intent(in) :: qinc
integer,intent(in) :: it ! guess time level
real(r_kind),dimension(:,: ),pointer:: ges_q2=>NULL()
integer(i_kind) i,j,ier
!*******************************************************************************
!
! 2-m temperature
! do it=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(it),'q2m',ges_q2,ier)
if(ier/=0) return
do j=1,lon2
do i=1,lat2
ges_q2(i,j) = ges_q2(i,j) + qinc(i,j)
end do
end do
! end do
return
end subroutine gsd_update_q2
subroutine gsd_gen_coast_prox
!$$$ subprogram documentation block
! . . . .
! subprogram: gsd_gen_coast_prox calculate coast proximity based on isli
! prgmmr: Hu org: GSD date: 2015-01-14
!
! abstract: This routine does the following things:
! 1) calculate coast proximity based on isli
!
!
! program history log:
! 2014-01-22 Hu - original code
!
! input argument list:
!
! output argument list:
!
! comments:
!
! attributes:
!$$$
use kinds, only: r_kind,i_kind
use mpimod, only: mype
use gridmod, only: lat2,lon2
use general_sub2grid_mod, only: general_gather2grid,general_scatter2sub
use general_commvars_mod, only: g1
use guess_grids, only: isli,coast_prox
use rapidrefresh_cldsurf_mod, only: l_gsd_soiltq_nudge
implicit none
! Declare passed variables
real(r_kind),dimension(:),allocatable:: worksub
real(r_kind),allocatable,dimension(:,:,:):: hwork
real(r_kind),allocatable,dimension(:,:,:):: hcoast_prox
integer(i_kind) workpe,ii
integer(i_kind) i,j,ico
integer(i_kind) ia,ib,ja,jb,ic,jc,nco,nip
!*******************************************************************************
!
if( l_gsd_soiltq_nudge) then
! water, land, seaice index
allocate(worksub(g1%inner_vars*g1%nlat*g1%nlon))
allocate(hwork(g1%inner_vars,g1%nlat,g1%nlon))
allocate(hcoast_prox(g1%inner_vars,g1%nlat,g1%nlon))
workpe=0
ii=0
do j=1,lon2
do i=1,lat2
ii=ii+1
worksub(ii)=isli(i,j,1)
end do
end do
call general_gather2grid(g1,worksub,hwork,workpe)
if(mype==workpe) then
!
! isli = 0 water, =1 land, =2 sea ice (on water)
hcoast_prox=0.0_r_kind
ico = 3
do j=1,g1%nlon
ja = max(1 ,j-ico)
jb = min(g1%nlon,j+ico+1)
do i=1,g1%nlat
if (abs(hwork(1,i,j)-1.0_r_kind) <0.001_r_kind .or. &
abs(hwork(1,i,j)-2.0_r_kind) <0.001_r_kind ) then
ia = max(1 ,i-ico)
ib = min(g1%nlat,i+ico+1)
nco = 0
nip = 0
do jc=ja,jb
do ic=ia,ib
if (abs(hwork(1,ic,jc)-1.0_r_kind) <0.001_r_kind .or. &
abs(hwork(1,ic,jc)-2.0_r_kind) <0.001_r_kind ) nco = nco+1
nip = nip+1
end do
end do
hcoast_prox(1,i,j) = float(nco)/float (nip)
end if
end do
end do
endif ! mype==workpe
deallocate(hwork)
!
worksub=0.0
call general_scatter2sub(g1,hcoast_prox,worksub,workpe)
deallocate(hcoast_prox)
ii=0
do j=1,lon2
do i=1,lat2
ii=ii+1
coast_prox(i,j)=worksub(ii)
end do
end do
else
coast_prox=0.0_r_kind
endif
return
end subroutine gsd_gen_coast_prox
end module gsd_update_mod