SCALE-RM
scale_cpl_phy_sfc_skin.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
21  !-----------------------------------------------------------------------------
22  implicit none
23  private
24  !-----------------------------------------------------------------------------
25  !
26  !++ Public procedure
27  !
28  public :: cpl_phy_sfc_skin_setup
29  public :: cpl_phy_sfc_skin
30 
31  !-----------------------------------------------------------------------------
32  !
33  !++ Public parameters & variables
34  !
35  !-----------------------------------------------------------------------------
36  !
37  !++ Private procedure
38  !
39  !-----------------------------------------------------------------------------
40  !
41  !++ Private parameters & variables
42  !
43  integer, private :: cpl_phy_sfc_skin_itr_max = 100 ! maximum iteration number
44 
45  real(RP), private :: cpl_phy_sfc_skin_dts_max = 5.0e-2_rp ! maximum delta surface temperature [K/s]
46  real(RP), private :: cpl_phy_sfc_skin_res_min = 1.0e+0_rp ! minimum value of residual
47  real(RP), private :: cpl_phy_sfc_skin_err_min = 1.0e-2_rp ! minimum value of error
48  real(RP), private :: cpl_phy_sfc_skin_dreslim = 1.0e+2_rp ! limiter of d(residual)
49 
50  logical, private :: initialized = .false.
51 
52  !-----------------------------------------------------------------------------
53 contains
54  !-----------------------------------------------------------------------------
56  subroutine cpl_phy_sfc_skin_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
90  end subroutine cpl_phy_sfc_skin_setup
91 
92  !-----------------------------------------------------------------------------
93  subroutine cpl_phy_sfc_skin( &
94  IA, IS, IE, &
95  JA, JS, JE, &
96  TMPA, PRSA, &
97  WA, UA, VA, &
98  RHOA, QVA, LH, &
99  Z1, PBL, &
100  RHOS, PRSS, &
101  RFLXD, &
102  TG, QVEF, &
103  ALBEDO, &
104  Rb, TC_dZ, &
105  Z0M, Z0H, Z0E, &
106  calc_flag, dt, &
107  model_name, &
108  TMPS, &
109  ZMFLX, XMFLX, YMFLX, &
110  SHFLX, QVFLX, GFLX, &
111  U10, V10, T2, Q2 )
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
499  end subroutine cpl_phy_sfc_skin
500 
501 end module scale_cpl_phy_sfc_skin
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
module atmosphere / saturation
real(rp), parameter, public const_stb
Stefan-Boltzman constant [W/m2/K4].
Definition: scale_const.F90:49
module coupler / surface-atmospehre
integer, parameter, public i_r_vis
module coupler / physics / surface skin
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
integer, parameter, public n_rad_dir
integer, parameter, public n_rad_rgn
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
module PROCESS
Definition: scale_prc.F90:11
procedure(bc), pointer, public bulkflux
subroutine, public cpl_phy_sfc_skin_setup
Setup.
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:89
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
module CONSTANT
Definition: scale_const.F90:11
integer, parameter, public i_r_direct
module profiler
Definition: scale_prof.F90:11
module Surface bulk flux
integer, parameter, public i_r_nir
module PRECISION
integer, parameter, public i_r_ir
module STDIO
Definition: scale_io.F90:10
integer, parameter, public i_r_diffuse
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)