SCALE-RM
Functions/Subroutines
scale_cpl_phy_sfc_skin Module Reference

module coupler / physics / surface skin More...

Functions/Subroutines

subroutine, public cpl_phy_sfc_skin_setup
 Setup. More...
 
subroutine, public cpl_phy_sfc_skin (IA, IS, IE, JA, JS, JE, TMPA, PRSA, WA, UA, VA, RHOA, QVA, LH, Z1, PBL, RHOS, PRSS, RFLXD, TG, QVEF, ALBEDO, Rb, TC_dZ, Z0M, Z0H, Z0E, calc_flag, dt, model_name, TMPS, ZMFLX, XMFLX, YMFLX, SHFLX, QVFLX, GFLX, U10, V10, T2, Q2)
 

Detailed Description

module coupler / physics / surface skin

Description
Skin surface model
Author
Team SCALE
NAMELIST
  • PARAM_CPL_PHY_SFC_SKIN
    nametypedefault valuecomment
    CPL_PHY_SFC_SKIN_ITR_MAX integer 100 maximum iteration number
    CPL_PHY_SFC_SKIN_DTS_MAX real(RP) 5.0E-2_RP maximum delta surface temperature [K/s]
    CPL_PHY_SFC_SKIN_RES_MIN real(RP) 1.0E+0_RP minimum value of residual
    CPL_PHY_SFC_SKIN_ERR_MIN real(RP) 1.0E-2_RP minimum value of error
    CPL_PHY_SFC_SKIN_DRESLIM real(RP) 1.0E+2_RP limiter of d(residual)

History Output
No history output

Function/Subroutine Documentation

◆ cpl_phy_sfc_skin_setup()

subroutine, public scale_cpl_phy_sfc_skin::cpl_phy_sfc_skin_setup ( )

Setup.

Definition at line 57 of file scale_cpl_phy_sfc_skin.F90.

References scale_io::io_fid_conf, and scale_prc::prc_abort().

Referenced by mod_land_driver::land_driver_setup().

57  use scale_prc, only: &
58  prc_abort
59  implicit none
60 
61  namelist / param_cpl_phy_sfc_skin / &
62  cpl_phy_sfc_skin_itr_max, &
63  cpl_phy_sfc_skin_dts_max, &
64  cpl_phy_sfc_skin_res_min, &
65  cpl_phy_sfc_skin_err_min, &
66  cpl_phy_sfc_skin_dreslim
67 
68  integer :: ierr
69  !---------------------------------------------------------------------------
70 
71  if ( initialized ) return
72 
73  log_newline
74  log_info("CPL_PHY_SFC_SKIN_setup",*) 'Setup'
75 
76  !--- read namelist
77  rewind(io_fid_conf)
78  read(io_fid_conf,nml=param_cpl_phy_sfc_skin,iostat=ierr)
79  if( ierr < 0 ) then !--- missing
80  log_info("CPL_PHY_SFC_SKIN_setup",*) 'Not found namelist. Default used.'
81  elseif( ierr > 0 ) then !--- fatal error
82  log_error("CPL_PHY_SFC_SKIN_setup",*) 'Not appropriate names in namelist PARAM_CPL_PHY_SFC_SKIN. Check!'
83  call prc_abort
84  endif
85  log_nml(param_cpl_phy_sfc_skin)
86 
87  initialized = .true.
88 
89  return
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cpl_phy_sfc_skin()

subroutine, public scale_cpl_phy_sfc_skin::cpl_phy_sfc_skin ( integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension (ia,ja), intent(in)  TMPA,
real(rp), dimension (ia,ja), intent(in)  PRSA,
real(rp), dimension (ia,ja), intent(in)  WA,
real(rp), dimension (ia,ja), intent(in)  UA,
real(rp), dimension (ia,ja), intent(in)  VA,
real(rp), dimension (ia,ja), intent(in)  RHOA,
real(rp), dimension (ia,ja), intent(in)  QVA,
real(rp), dimension (ia,ja), intent(in)  LH,
real(rp), dimension (ia,ja), intent(in)  Z1,
real(rp), dimension (ia,ja), intent(in)  PBL,
real(rp), dimension (ia,ja), intent(in)  RHOS,
real(rp), dimension (ia,ja), intent(in)  PRSS,
real(rp), dimension (ia,ja,n_rad_dir,n_rad_rgn), intent(in)  RFLXD,
real(rp), dimension (ia,ja), intent(in)  TG,
real(rp), dimension (ia,ja), intent(in)  QVEF,
real(rp), dimension (ia,ja,n_rad_dir,n_rad_rgn), intent(in)  ALBEDO,
real(rp), dimension (ia,ja), intent(in)  Rb,
real(rp), dimension (ia,ja), intent(in)  TC_dZ,
real(rp), dimension (ia,ja), intent(in)  Z0M,
real(rp), dimension (ia,ja), intent(in)  Z0H,
real(rp), dimension (ia,ja), intent(in)  Z0E,
logical, dimension(ia,ja), intent(in)  calc_flag,
real(dp), intent(in)  dt,
character(len=*), intent(in)  model_name,
real(rp), dimension (ia,ja), intent(inout)  TMPS,
real(rp), dimension (ia,ja), intent(out)  ZMFLX,
real(rp), dimension (ia,ja), intent(out)  XMFLX,
real(rp), dimension (ia,ja), intent(out)  YMFLX,
real(rp), dimension (ia,ja), intent(out)  SHFLX,
real(rp), dimension (ia,ja), intent(out)  QVFLX,
real(rp), dimension (ia,ja), intent(out)  GFLX,
real(rp), dimension (ia,ja), intent(out)  U10,
real(rp), dimension (ia,ja), intent(out)  V10,
real(rp), dimension (ia,ja), intent(out)  T2,
real(rp), dimension (ia,ja), intent(out)  Q2 
)

Definition at line 112 of file scale_cpl_phy_sfc_skin.F90.

References scale_bulkflux::bulkflux, scale_const::const_cpdry, scale_const::const_pre00, scale_const::const_rdry, scale_const::const_rvap, scale_const::const_stb, scale_cpl_sfc_index::i_r_diffuse, scale_cpl_sfc_index::i_r_direct, scale_cpl_sfc_index::i_r_ir, scale_cpl_sfc_index::i_r_nir, scale_cpl_sfc_index::i_r_vis, scale_prc::prc_abort(), and scale_prc::prc_myrank.

Referenced by mod_land_driver::land_driver_calc_tendency().

112  use scale_prc, only: &
113  prc_myrank, &
114  prc_abort
115  use scale_const, only: &
116  pre00 => const_pre00, &
117  rdry => const_rdry, &
118  cpdry => const_cpdry, &
119  rvap => const_rvap, &
120  stb => const_stb
121  use scale_atmos_saturation, only: &
122  qsat => atmos_saturation_pres2qsat_all
123  use scale_bulkflux, only: &
124  bulkflux
125  implicit none
126 
127  integer, intent(in) :: ia, is, ie
128  integer, intent(in) :: ja, js, je
129  real(RP), intent(in) :: tmpa (ia,ja) ! temperature at the lowest atmospheric layer [K]
130  real(RP), intent(in) :: prsa (ia,ja) ! pressure at the lowest atmospheric layer [Pa]
131  real(RP), intent(in) :: wa (ia,ja) ! velocity w at the lowest atmospheric layer [m/s]
132  real(RP), intent(in) :: ua (ia,ja) ! velocity u at the lowest atmospheric layer [m/s]
133  real(RP), intent(in) :: va (ia,ja) ! velocity v at the lowest atmospheric layer [m/s]
134  real(RP), intent(in) :: rhoa (ia,ja) ! density at the lowest atmospheric layer [kg/m3]
135  real(RP), intent(in) :: qva (ia,ja) ! ratio of water vapor mass to total mass at the lowest atmospheric layer [kg/kg]
136  real(RP), intent(in) :: lh (ia,ja) ! latent heat at the lowest atmospheric layer [J/kg]
137  real(RP), intent(in) :: z1 (ia,ja) ! cell center height at the lowest atmospheric layer [m]
138  real(RP), intent(in) :: pbl (ia,ja) ! the top of atmospheric mixing layer [m]
139  real(RP), intent(in) :: rhos (ia,ja) ! density at the surface [kg/m3]
140  real(RP), intent(in) :: prss (ia,ja) ! pressure at the surface [Pa]
141  real(RP), intent(in) :: rflxd (ia,ja,n_rad_dir,n_rad_rgn) ! downward radiation flux at the surface (direct/diffuse,IR/near-IR/VIS) [J/m2/s]
142  real(RP), intent(in) :: tg (ia,ja) ! subsurface temperature [K]
143  real(RP), intent(in) :: qvef (ia,ja) ! efficiency of evaporation (0-1)
144  real(RP), intent(in) :: albedo (ia,ja,n_rad_dir,n_rad_rgn) ! surface albedo (direct/diffuse,IR/near-IR/VIS) (0-1)
145  real(RP), intent(in) :: rb (ia,ja) ! stomata resistance [1/s]
146  real(RP), intent(in) :: tc_dz (ia,ja) ! thermal conductivity / depth between surface and subsurface [J/m2/s/K]
147  real(RP), intent(in) :: z0m (ia,ja) ! roughness length for momemtum [m]
148  real(RP), intent(in) :: z0h (ia,ja) ! roughness length for heat [m]
149  real(RP), intent(in) :: z0e (ia,ja) ! roughness length for vapor [m]
150  logical, intent(in) :: calc_flag(ia,ja) ! to decide calculate or not
151  real(DP), intent(in) :: dt ! delta time
152  character(len=*), intent(in) :: model_name
153  real(RP), intent(inout) :: tmps (ia,ja) ! surface temperature [K]
154  real(RP), intent(out) :: zmflx (ia,ja) ! z-momentum flux at the surface [kg/m/s2]
155  real(RP), intent(out) :: xmflx (ia,ja) ! x-momentum flux at the surface [kg/m/s2]
156  real(RP), intent(out) :: ymflx (ia,ja) ! y-momentum flux at the surface [kg/m/s2]
157  real(RP), intent(out) :: shflx (ia,ja) ! sensible heat flux at the surface [J/m2/s]
158  real(RP), intent(out) :: qvflx (ia,ja) ! water vapor flux at the surface [kg/m2/s]
159  real(RP), intent(out) :: gflx (ia,ja) ! subsurface heat flux at the surface [J/m2/s]
160  real(RP), intent(out) :: u10 (ia,ja) ! velocity u at 10m [m/s]
161  real(RP), intent(out) :: v10 (ia,ja) ! velocity v at 10m [m/s]
162  real(RP), intent(out) :: t2 (ia,ja) ! temperature at 2m [K]
163  real(RP), intent(out) :: q2 (ia,ja) ! water vapor at 2m [kg/kg]
164 
165  real(RP), parameter :: dts0 = 1.0e-4_rp ! delta surface temp.
166  real(RP), parameter :: redf_min = 1.0e-2_rp ! minimum reduced factor
167  real(RP), parameter :: redf_max = 1.0e+0_rp ! maximum reduced factor
168  real(RP), parameter :: tfa = 0.5e+0_rp ! factor a in Tomita (2009)
169  real(RP), parameter :: tfb = 1.1e+0_rp ! factor b in Tomita (2009)
170 
171  real(RP) :: tmps1(ia,ja)
172 
173  real(RP) :: emis ! surface longwave emission [J/m2/s]
174  real(RP) :: lwd ! surface downward longwave radiation flux [J/m2/s]
175  real(RP) :: lwu ! surface upward longwave radiation flux [J/m2/s]
176  real(RP) :: swd ! surface downward shortwave radiation flux [J/m2/s]
177  real(RP) :: swu ! surface upward shortwave radiation flux [J/m2/s]
178  real(RP) :: res ! residual
179 
180  real(RP) :: dres ! d(residual)/dTMPS
181  real(RP) :: oldres ! residual in previous step
182  real(RP) :: redf ! reduced factor
183 
184  real(RP) :: ustar, dustar ! friction velocity [m]
185  real(RP) :: tstar, dtstar ! friction potential temperature [K]
186  real(RP) :: qstar, dqstar ! friction water vapor mass ratio [kg/kg]
187  real(RP) :: uabs, duabs ! modified absolute velocity [m/s]
188  real(RP) :: ra, dra ! Aerodynamic resistance (=1/Ce) [1/s]
189 
190  real(RP) :: qvsat, dqvsat ! saturation water vapor mixing ratio at surface [kg/kg]
191  real(RP) :: qvs, dqvs ! water vapor mixing ratio at surface [kg/kg]
192  real(RP) :: rtot ! total gas constant
193  real(RP) :: qdry ! dry air mass ratio [kg/kg]
194 
195  real(RP) :: fracu10 ! calculation parameter for U10 [1]
196  real(RP) :: fract2 ! calculation parameter for T2 [1]
197  real(RP) :: fracq2 ! calculation parameter for Q2 [1]
198 
199  integer :: i, j, n
200  !---------------------------------------------------------------------------
201 
202  log_progress(*) 'coupler / physics / surface / SKIN'
203 
204  ! copy surfce temperature for iteration
205  !$omp parallel do
206  do j = js, je
207  do i = is, ie
208  tmps1(i,j) = tmps(i,j)
209  enddo
210  enddo
211 
212  ! update surface temperature
213 #ifndef __GFORTRAN__
214  !$omp parallel do default(none) &
215  !$omp private(qdry,Rtot,redf,res,emis,LWD,LWU,SWD,SWU,dres,oldres,QVS,dQVS, &
216  !$omp QVsat,dQVsat,Ustar,dUstar,Tstar,dTstar,Qstar,dQstar,Uabs,dUabs,Ra,dRa,FracU10,FracT2,FracQ2) &
217  !$omp shared(IS,IE,JS,JE,Rdry,CPdry,PRC_myrank,IO_FID_LOG,IO_L,model_name,bulkflux, &
218  !$omp CPL_PHY_SFC_SKIN_itr_max,CPL_PHY_SFC_SKIN_dTS_max,CPL_PHY_SFC_SKIN_dreslim,CPL_PHY_SFC_SKIN_err_min, CPL_PHY_SFC_SKIN_res_min, &
219  !$omp calc_flag,dt,QVA,TMPA,PRSA,RHOA,WA,UA,VA,LH,Z1,PBL, &
220  !$omp TG,PRSS,RHOS,TMPS1,QVEF,Z0M,Z0H,Z0E,Rb,TC_dZ,ALBEDO,RFLXD, &
221  !$omp TMPS,ZMFLX,XMFLX,YMFLX,SHFLX,QVFLX,GFLX,U10,V10,T2,Q2)
222 #else
223  !$omp parallel do default(shared) &
224  !$omp private(qdry,Rtot,redf,res,emis,LWD,LWU,SWD,SWU,dres,oldres,QVS,dQVS, &
225  !$omp QVsat,dQVsat,Ustar,dUstar,Tstar,dTstar,Qstar,dQstar,Uabs,dUabs,Ra,dRa,FracU10,FracT2,FracQ2)
226 #endif
227  do j = js, je
228  do i = is, ie
229  if ( calc_flag(i,j) ) then
230 
231  qdry = 1.0_rp - qva(i,j)
232  rtot = qdry * rdry + qva(i,j) * rvap
233 
234  redf = 1.0_rp
235  oldres = huge(0.0_rp)
236 
237  ! modified Newton-Raphson method (Tomita 2009)
238  do n = 1, cpl_phy_sfc_skin_itr_max
239 
240  call qsat( tmps1(i,j), prss(i,j), qdry, qvsat )
241  call qsat( tmps1(i,j)+dts0, prss(i,j), qdry, dqvsat )
242 
243  qvs = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
244  + ( qvef(i,j) ) * qvsat
245  dqvs = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
246  + ( qvef(i,j) ) * dqvsat
247 
248  call bulkflux( ustar, & ! [OUT]
249  tstar, & ! [OUT]
250  qstar, & ! [OUT]
251  uabs, & ! [OUT]
252  ra, & ! [OUT]
253  fracu10, & ! [OUT] ! not used
254  fract2, & ! [OUT] ! not used
255  fracq2, & ! [OUT] ! not used
256  tmpa(i,j), & ! [IN]
257  tmps1(i,j), & ! [IN]
258  prsa(i,j), & ! [IN]
259  prss(i,j), & ! [IN]
260  qva(i,j), & ! [IN]
261  qvs, & ! [IN]
262  ua(i,j), & ! [IN]
263  va(i,j), & ! [IN]
264  z1(i,j), & ! [IN]
265  pbl(i,j), & ! [IN]
266  z0m(i,j), & ! [IN]
267  z0h(i,j), & ! [IN]
268  z0e(i,j) ) ! [IN]
269 
270  call bulkflux( dustar, & ! [OUT]
271  dtstar, & ! [OUT]
272  dqstar, & ! [OUT]
273  duabs, & ! [OUT]
274  dra, & ! [OUT] ! not used
275  fracu10, & ! [OUT] ! not used
276  fract2, & ! [OUT] ! not used
277  fracq2, & ! [OUT] ! not used
278  tmpa(i,j), & ! [IN]
279  tmps1(i,j)+dts0, & ! [IN]
280  prsa(i,j), & ! [IN]
281  prss(i,j), & ! [IN]
282  qva(i,j), & ! [IN]
283  dqvs, & ! [IN]
284  ua(i,j), & ! [IN]
285  va(i,j), & ! [IN]
286  z1(i,j), & ! [IN]
287  pbl(i,j), & ! [IN]
288  z0m(i,j), & ! [IN]
289  z0h(i,j), & ! [IN]
290  z0e(i,j) ) ! [IN]
291 
292  emis = ( 1.0_rp - albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps1(i,j)**4
293 
294  lwd = rflxd(i,j,i_r_diffuse,i_r_ir)
295  lwu = rflxd(i,j,i_r_diffuse,i_r_ir) * albedo(i,j,i_r_diffuse,i_r_ir) + emis
296  swd = rflxd(i,j,i_r_direct ,i_r_nir) &
297  + rflxd(i,j,i_r_diffuse,i_r_nir) &
298  + rflxd(i,j,i_r_direct ,i_r_vis) &
299  + rflxd(i,j,i_r_diffuse,i_r_vis)
300  swu = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
301  + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
302  + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
303  + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)
304 
305  ! calculation for residual
306  res = swd - swu + lwd - lwu &
307  + cpdry * rhos(i,j) * ustar * tstar &
308  + lh(i,j) * rhos(i,j) * ustar * qstar * ra / ( ra+rb(i,j) ) &
309  - tc_dz(i,j) * ( tmps1(i,j) - tg(i,j) )
310 
311  ! calculation for d(residual)/dTMPS
312  dres = -4.0_rp * emis / tmps1(i,j) &
313  + cpdry * rhos(i,j) * ( ustar*(dtstar-tstar)/dts0 + tstar*(dustar-ustar)/dts0 ) &
314  + lh(i,j) * rhos(i,j) * ( ustar*(dqstar-qstar)/dts0 + qstar*(dustar-ustar)/dts0 ) * ra / ( ra+rb(i,j) ) &
315  - tc_dz(i,j)
316 
317  ! convergence test with residual and error levels
318  if ( abs(res ) < cpl_phy_sfc_skin_res_min &
319  .OR. abs(res/dres) < cpl_phy_sfc_skin_err_min ) then
320  exit
321  endif
322 
323  ! stop iteration to prevent numerical error
324  if ( abs(dres) * cpl_phy_sfc_skin_dreslim < abs(res) ) then
325  exit
326  endif
327 
328  ! calculate reduced factor
329  if ( dres < 0.0_rp ) then
330  if ( abs(res) > abs(oldres) ) then
331  redf = max( tfa*abs(redf), redf_min )
332  else
333  redf = min( tfb*abs(redf), redf_max )
334  endif
335  else
336  redf = -1.0_rp
337  endif
338 
339  ! estimate next surface temperature
340  tmps1(i,j) = tmps1(i,j) - redf * res / dres
341 
342  ! save residual in this step
343  oldres = res
344  enddo
345 
346  ! update surface temperature with limitation
347  tmps1(i,j) = min( max( tmps1(i,j), &
348  tmps(i,j) - cpl_phy_sfc_skin_dts_max * real(dt,kind=RP) ), &
349  tmps (i,j) + cpl_phy_sfc_skin_dts_max * real(dt,kind=RP) )
350 
351  if ( n > cpl_phy_sfc_skin_itr_max ) then
352  ! surface temperature was not converged
353  log_warn("CPL_PHY_SFC_skin",*) 'surface tempearture was not converged. ', trim(model_name)
354  log_newline
355  log_info_cont('(A,I32)' ) 'PRC_myrank [no unit] :', prc_myrank
356  log_info_cont('(A,I32)' ) 'number of i [no unit] :', i
357  log_info_cont('(A,I32)' ) 'number of j [no unit] :', j
358  log_newline
359  log_info_cont('(A,I32)' ) 'loop number [no unit] :', n
360  log_info_cont('(A,F32.16)') 'Residual [J/m2/s] :', res
361  log_info_cont('(A,F32.16)') 'delta Residual [J/m2/s] :', dres
362  log_newline
363  log_info_cont('(A,F32.16)') 'temperature [K] :', tmpa(i,j)
364  log_info_cont('(A,F32.16)') 'pressure [Pa] :', prsa(i,j)
365  log_info_cont('(A,F32.16)') 'velocity w [m/s] :', wa(i,j)
366  log_info_cont('(A,F32.16)') 'velocity u [m/s] :', ua(i,j)
367  log_info_cont('(A,F32.16)') 'velocity v [m/s] :', va(i,j)
368  log_info_cont('(A,F32.16)') 'density [kg/m3] :', rhoa(i,j)
369  log_info_cont('(A,F32.16)') 'water vapor mass ratio [kg/kg] :', qva(i,j)
370  log_info_cont('(A,F32.16)') 'cell center height [m] :', z1(i,j)
371  log_info_cont('(A,F32.16)') 'atmospheric mixing layer height [m] :', pbl(i,j)
372  log_info_cont('(A,F32.16)') 'pressure at the surface [Pa] :', prss(i,j)
373  log_info_cont('(A,F32.16)') 'downward radiation (IR, direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_ir )
374  log_info_cont('(A,F32.16)') 'downward radiation (IR, diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_ir )
375  log_info_cont('(A,F32.16)') 'downward radiation (NIR,direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_nir)
376  log_info_cont('(A,F32.16)') 'downward radiation (NIR,diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_nir)
377  log_info_cont('(A,F32.16)') 'downward radiation (VIS,direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_vis)
378  log_info_cont('(A,F32.16)') 'downward radiation (VIS,diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_vis)
379  log_newline
380  log_info_cont('(A,F32.16)') 'soil temperature [K] :', tg(i,j)
381  log_info_cont('(A,F32.16)') 'surface temperature [K] :', tmps(i,j)
382  log_info_cont('(A,F32.16)') 'efficiency of evaporation [1] :', qvef(i,j)
383  log_info_cont('(A,F32.16)') 'surface albedo (IR, direct ) [1] :', albedo(i,j,i_r_direct ,i_r_ir )
384  log_info_cont('(A,F32.16)') 'surface albedo (IR, diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_ir )
385  log_info_cont('(A,F32.16)') 'surface albedo (NIR,direct ) [1] :', albedo(i,j,i_r_direct ,i_r_nir)
386  log_info_cont('(A,F32.16)') 'surface albedo (NIR,diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_nir)
387  log_info_cont('(A,F32.16)') 'surface albedo (VIS,direct ) [1] :', albedo(i,j,i_r_direct ,i_r_vis)
388  log_info_cont('(A,F32.16)') 'surface albedo (VIS,diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_vis)
389  log_info_cont('(A,F32.16)') 'thermal conductivity / depth [J/m2/s/K] :', tc_dz(i,j)
390  log_info_cont('(A,F32.16)') 'roughness length for momemtum [m] :', z0m(i,j)
391  log_info_cont('(A,F32.16)') 'roughness length for heat [m] :', z0h(i,j)
392  log_info_cont('(A,F32.16)') 'roughness length for vapor [m] :', z0e(i,j)
393  log_newline
394  log_info_cont('(A,F32.16)') 'latent heat [J/kg] :', lh(i,j)
395  log_info_cont('(A,F32.16)') 'friction velocity [m] :', ustar
396  log_info_cont('(A,F32.16)') 'friction potential temperature [K] :', tstar
397  log_info_cont('(A,F32.16)') 'friction water vapor mass ratio [kg/kg] :', qstar
398  log_info_cont('(A,F32.16)') 'd(friction velocity) [m] :', dustar
399  log_info_cont('(A,F32.16)') 'd(friction potential temperature) [K] :', dtstar
400  log_info_cont('(A,F32.16)') 'd(friction water vapor mass ratio) [kg/kg] :', dqstar
401  log_info_cont('(A,F32.16)') 'modified absolute velocity [m/s] :', uabs
402  log_info_cont('(A,F32.16)') 'next surface temperature [K] :', tmps1(i,j)
403 
404  ! check NaN
405  if ( .NOT. ( res > -1.0_rp .OR. res < 1.0_rp ) ) then ! must be NaN
406  log_error("CPL_PHY_SFC_skin",*) 'NaN is detected for surface temperature. ', trim(model_name)
407  call prc_abort
408  endif
409  endif
410 
411  ! calculate surface flux
412  tmps(i,j) = tmps1(i,j)
413 
414  qdry = 1.0_rp - qva(i,j)
415  rtot = qdry * rdry + qva(i,j) * rvap
416 
417  call qsat( tmps(i,j), prss(i,j), qdry, qvsat )
418 
419  qvs = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
420  + ( qvef(i,j) ) * qvsat
421 
422  call bulkflux( ustar, & ! [OUT]
423  tstar, & ! [OUT]
424  qstar, & ! [OUT]
425  uabs, & ! [OUT]
426  ra, & ! [OUT]
427  fracu10, & ! [OUT]
428  fract2, & ! [OUT]
429  fracq2, & ! [OUT]
430  tmpa(i,j), & ! [IN]
431  tmps(i,j), & ! [IN]
432  prsa(i,j), & ! [IN]
433  prss(i,j), & ! [IN]
434  qva(i,j), & ! [IN]
435  qvs, & ! [IN]
436  ua(i,j), & ! [IN]
437  va(i,j), & ! [IN]
438  z1(i,j), & ! [IN]
439  pbl(i,j), & ! [IN]
440  z0m(i,j), & ! [IN]
441  z0h(i,j), & ! [IN]
442  z0e(i,j) ) ! [IN]
443 
444  zmflx(i,j) = -rhos(i,j) * ustar * ustar / uabs * wa(i,j)
445  xmflx(i,j) = -rhos(i,j) * ustar * ustar / uabs * ua(i,j)
446  ymflx(i,j) = -rhos(i,j) * ustar * ustar / uabs * va(i,j)
447  shflx(i,j) = -rhos(i,j) * ustar * tstar * cpdry
448  qvflx(i,j) = -rhos(i,j) * ustar * qstar * ra / ( ra+rb(i,j) )
449 
450  emis = ( 1.0_rp-albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps(i,j)**4
451 
452  lwd = rflxd(i,j,i_r_diffuse,i_r_ir)
453  lwu = rflxd(i,j,i_r_diffuse,i_r_ir) * albedo(i,j,i_r_diffuse,i_r_ir) + emis
454  swd = rflxd(i,j,i_r_direct ,i_r_nir) &
455  + rflxd(i,j,i_r_diffuse,i_r_nir) &
456  + rflxd(i,j,i_r_direct ,i_r_vis) &
457  + rflxd(i,j,i_r_diffuse,i_r_vis)
458  swu = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
459  + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
460  + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
461  + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)
462 
463  gflx(i,j) = -tc_dz(i,j) * ( tmps(i,j) - tg(i,j) )
464 
465  ! calculation for residual
466  res = swd - swu + lwd - lwu - shflx(i,j) - qvflx(i,j) * lh(i,j) + gflx(i,j)
467 
468  ! put residual in ground heat flux
469  gflx(i,j) = gflx(i,j) - res
470 
471  ! diagnostic variables considering unstable/stable state
472  !U10(i,j) = FracU10 * UA(i,j)
473  !V10(i,j) = FracU10 * VA(i,j)
474  !T2 (i,j) = ( 1.0_RP - FracT2 ) * TMPS(i,j) + FracT2 * TMPA(i,j)
475  !Q2 (i,j) = ( 1.0_RP - FracQ2 ) * QVS + FracQ2 * QVA (i,j)
476 
477  ! diagnostic variables for neutral state
478  u10(i,j) = ua(i,j) * log( 10.0_rp / z0m(i,j) ) / log( z1(i,j) / z0m(i,j) )
479  v10(i,j) = va(i,j) * log( 10.0_rp / z0m(i,j) ) / log( z1(i,j) / z0m(i,j) )
480  t2(i,j) = tmps(i,j) + ( tmpa(i,j) - tmps(i,j) ) * log( 2.0_rp / z0h(i,j) ) / log( z1(i,j) / z0h(i,j) )
481  q2(i,j) = qvs + ( qva(i,j) - qvs ) * log( 2.0_rp / z0e(i,j) ) / log( z1(i,j) / z0e(i,j) )
482 
483  else ! not calculate surface flux
484  zmflx(i,j) = 0.0_rp
485  xmflx(i,j) = 0.0_rp
486  ymflx(i,j) = 0.0_rp
487  shflx(i,j) = 0.0_rp
488  qvflx(i,j) = 0.0_rp
489  gflx(i,j) = 0.0_rp
490  u10(i,j) = 0.0_rp
491  v10(i,j) = 0.0_rp
492  t2(i,j) = 0.0_rp
493  q2(i,j) = 0.0_rp
494  endif
495  enddo
496  enddo
497 
498  return
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_index::va
integer, public va
Definition: scale_index.F90:35
scale_const::const_stb
real(rp), parameter, public const_stb
Stefan-Boltzman constant [W/m2/K4].
Definition: scale_const.F90:49
scale_atmos_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_grid_cartesc_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_bulkflux::bulkflux
procedure(bc), pointer, public bulkflux
Definition: scale_bulkflux.F90:80
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:89
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_bulkflux
module Surface bulk flux
Definition: scale_bulkflux.F90:12
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by   doxygen 1.8.13