SCALE-RM
scale_cpl_phy_sfc_skin.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_cpl_phy_sfc_skin
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  use scale_cpl_sfc_index
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 
49  logical, private :: initialized = .false.
50 
51  !-----------------------------------------------------------------------------
52 contains
53  !-----------------------------------------------------------------------------
55  subroutine cpl_phy_sfc_skin_setup
56  use scale_prc, only: &
57  prc_abort
58  implicit none
59 
60  namelist / param_cpl_phy_sfc_skin / &
61  cpl_phy_sfc_skin_itr_max, &
62  cpl_phy_sfc_skin_dts_max, &
63  cpl_phy_sfc_skin_res_min, &
64  cpl_phy_sfc_skin_err_min
65 
66  integer :: ierr
67  !---------------------------------------------------------------------------
68 
69  if ( initialized ) return
70 
71  log_newline
72  log_info("CPL_PHY_SFC_SKIN_setup",*) 'Setup'
73 
74  !--- read namelist
75  rewind(io_fid_conf)
76  read(io_fid_conf,nml=param_cpl_phy_sfc_skin,iostat=ierr)
77  if( ierr < 0 ) then !--- missing
78  log_info("CPL_PHY_SFC_SKIN_setup",*) 'Not found namelist. Default used.'
79  elseif( ierr > 0 ) then !--- fatal error
80  log_error("CPL_PHY_SFC_SKIN_setup",*) 'Not appropriate names in namelist PARAM_CPL_PHY_SFC_SKIN. Check!'
81  call prc_abort
82  endif
83  log_nml(param_cpl_phy_sfc_skin)
84 
85  initialized = .true.
86 
87  return
88  end subroutine cpl_phy_sfc_skin_setup
89 
90  !-----------------------------------------------------------------------------
91  subroutine cpl_phy_sfc_skin( &
92  IA, IS, IE, &
93  JA, JS, JE, &
94  TMPA, PRSA, &
95  WA, UA, VA, &
96  RHOA, QVA, &
97  LH, Z1, PBL, &
98  RHOS, PRSS, &
99  RFLXD, &
100  TG, WSTR, QVEF, &
101  ALBEDO, &
102  Rb, TC_dZ, &
103  Z0M, Z0H, Z0E, &
104  calc_flag, dt, &
105  model_name, &
106  TMPS, &
107  ZMFLX, XMFLX, YMFLX, &
108  SHFLX, LHFLX, QVFLX, &
109  GFLX, &
110  Ustar, Tstar, Qstar, &
111  Wstar, &
112  RLmo, &
113  U10, V10, T2, Q2 )
114  use scale_prc, only: &
115  prc_myrank, &
116  prc_abort
117  use scale_const, only: &
118  eps => const_eps, &
119  undef => const_undef, &
120  pre00 => const_pre00, &
121  tem00 => const_tem00, &
122  rdry => const_rdry, &
123  cpdry => const_cpdry, &
124  rvap => const_rvap, &
125  stb => const_stb
126  use scale_atmos_saturation, only: &
127  qsat => atmos_saturation_dens2qsat_all
128 ! qsat => ATMOS_SATURATION_pres2qsat_all
129  use scale_atmos_hydrometeor, only: &
130  atmos_hydrometeor_lhv, &
131  atmos_hydrometeor_lhs, &
132  cv_water, &
133  cv_ice, &
134  lhf
135  use scale_bulkflux, only: &
136  bulkflux, &
137  bulkflux_diagnose_surface
138  implicit none
139 
140  integer, intent(in) :: ia, is, ie
141  integer, intent(in) :: ja, js, je
142  real(rp), intent(in) :: tmpa (ia,ja) ! temperature at the lowest atmospheric layer [K]
143  real(rp), intent(in) :: prsa (ia,ja) ! pressure at the lowest atmospheric layer [Pa]
144  real(rp), intent(in) :: wa (ia,ja) ! velocity w at the lowest atmospheric layer [m/s]
145  real(rp), intent(in) :: ua (ia,ja) ! velocity u at the lowest atmospheric layer [m/s]
146  real(rp), intent(in) :: va (ia,ja) ! velocity v at the lowest atmospheric layer [m/s]
147  real(rp), intent(in) :: rhoa (ia,ja) ! density at the lowest atmospheric layer [kg/m3]
148  real(rp), intent(in) :: qva (ia,ja) ! ratio of water vapor mass to total mass at the lowest atmospheric layer [kg/kg]
149  real(rp), intent(in) :: lh (ia,ja) ! latent heat [J/kg]
150  real(rp), intent(in) :: z1 (ia,ja) ! cell center height at the lowest atmospheric layer [m]
151  real(rp), intent(in) :: pbl (ia,ja) ! the top of atmospheric mixing layer [m]
152  real(rp), intent(in) :: rhos (ia,ja) ! density at the surface [kg/m3]
153  real(rp), intent(in) :: prss (ia,ja) ! pressure at the surface [Pa]
154  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]
155  real(rp), intent(in) :: tg (ia,ja) ! subsurface temperature [K]
156  real(rp), intent(in) :: wstr (ia,ja) ! amount of water storage [kg/m2]
157  real(rp), intent(in) :: qvef (ia,ja) ! efficiency of evaporation (0-1)
158  real(rp), intent(in) :: albedo (ia,ja,n_rad_dir,n_rad_rgn) ! surface albedo (direct/diffuse,IR/near-IR/VIS) (0-1)
159  real(rp), intent(in) :: rb (ia,ja) ! stomata resistance [1/s]
160  real(rp), intent(in) :: tc_dz (ia,ja) ! thermal conductivity / depth between surface and subsurface [J/m2/s/K]
161  real(rp), intent(in) :: z0m (ia,ja) ! roughness length for momemtum [m]
162  real(rp), intent(in) :: z0h (ia,ja) ! roughness length for heat [m]
163  real(rp), intent(in) :: z0e (ia,ja) ! roughness length for vapor [m]
164  logical, intent(in) :: calc_flag(ia,ja) ! to decide calculate or not
165  real(dp), intent(in) :: dt ! delta time
166  character(len=*), intent(in) :: model_name
167 
168  real(rp), intent(inout) :: tmps (ia,ja) ! surface temperature [K]
169 
170  real(rp), intent(out) :: zmflx (ia,ja) ! z-momentum flux at the surface [kg/m/s2]
171  real(rp), intent(out) :: xmflx (ia,ja) ! x-momentum flux at the surface [kg/m/s2]
172  real(rp), intent(out) :: ymflx (ia,ja) ! y-momentum flux at the surface [kg/m/s2]
173  real(rp), intent(out) :: shflx (ia,ja) ! sensible heat flux at the surface [J/m2/s]
174  real(rp), intent(out) :: lhflx (ia,ja) ! latent heat flux at the surface [J/m2/s]
175  real(rp), intent(out) :: qvflx (ia,ja) ! water vapor flux at the surface [kg/m2/s]
176  real(rp), intent(out) :: gflx (ia,ja) ! subsurface heat flux at the surface [J/m2/s]
177  real(rp), intent(out) :: ustar (ia,ja) ! friction velocity [m/s]
178  real(rp), intent(out) :: tstar (ia,ja) ! temperature scale [K]
179  real(rp), intent(out) :: qstar (ia,ja) ! moisture scale [kg/kg]
180  real(rp), intent(out) :: wstar (ia,ja) ! convective velocity scale [m/s]
181  real(rp), intent(out) :: rlmo (ia,ja) ! inversed Obukhov length [1/m]
182  real(rp), intent(out) :: u10 (ia,ja) ! velocity u at 10m [m/s]
183  real(rp), intent(out) :: v10 (ia,ja) ! velocity v at 10m [m/s]
184  real(rp), intent(out) :: t2 (ia,ja) ! temperature at 2m [K]
185  real(rp), intent(out) :: q2 (ia,ja) ! water vapor at 2m [kg/kg]
186 
187  real(rp), parameter :: dts0 = 1.0e-4_rp ! delta surface temp.
188  real(rp), parameter :: redf_min = 1.0e-2_rp ! minimum reduced factor
189  real(rp), parameter :: redf_max = 1.0e+0_rp ! maximum reduced factor
190  real(rp), parameter :: tfa = 0.5e+0_rp ! factor a in Tomita (2009)
191  real(rp), parameter :: tfb = 1.1e+0_rp ! factor b in Tomita (2009)
192 
193  real(rp) :: tmps1(ia,ja)
194 
195  real(rp) :: emis ! surface longwave emission [J/m2/s]
196  real(rp) :: lwd ! surface downward longwave radiation flux [J/m2/s]
197  real(rp) :: lwu ! surface upward longwave radiation flux [J/m2/s]
198  real(rp) :: swd ! surface downward shortwave radiation flux [J/m2/s]
199  real(rp) :: swu ! surface upward shortwave radiation flux [J/m2/s]
200  real(rp) :: flx_qv ! surface upward qv flux [kg/m2/s]
201  real(rp) :: res ! residual
202 
203  real(rp) :: dres ! d(residual)/dTMPS
204  real(rp) :: oldres ! residual in previous step
205  real(rp) :: redf ! reduced factor
206  real(rp) :: dts ! temperature change
207  real(rp) :: olddts ! temperature change in previous step
208 
209  real(rp) :: dustar ! friction velocity difference [m/s]
210  real(rp) :: dtstar ! friction potential temperature difference [K]
211  real(rp) :: dqstar ! friction water vapor mass ratio difference [kg/kg]
212  real(rp) :: dwstar ! free convection velocity scale difference [m/s]
213  real(rp) :: drlmo ! inversed Obukhov length [1/m]
214  real(rp) :: uabs, duabs ! modified absolute velocity [m/s]
215  real(rp) :: ra, dra ! Aerodynamic resistance (=1/Ce) [1/s]
216 
217  real(rp) :: qvsat, dqvsat ! saturation water vapor mixing ratio at surface [kg/kg]
218  real(rp) :: qvs(ia,ja), dqvs ! water vapor mixing ratio at surface [kg/kg]
219  real(rp) :: rtot ! total gas constant
220  real(rp) :: qdry ! dry air mass ratio [kg/kg]
221 
222  real(rp) :: fracu10(ia,ja), dfracu10 ! calculation parameter for U10 [1]
223  real(rp) :: fract2 (ia,ja), dfract2 ! calculation parameter for T2 [1]
224  real(rp) :: fracq2 (ia,ja), dfracq2 ! calculation parameter for Q2 [1]
225 
226  real(rp) :: mflux
227 
228  integer :: i, j, n
229  !---------------------------------------------------------------------------
230 
231  log_progress(*) 'coupler / physics / surface / SKIN'
232 
233  ! copy surfce temperature for iteration
234  !$omp parallel do
235  do j = js, je
236  do i = is, ie
237  tmps1(i,j) = tmps(i,j)
238  enddo
239  enddo
240 
241  ! update surface temperature
242  !$omp parallel do &
243 #ifndef __GFORTRAN__
244  !$omp default(none) &
245  !$omp shared(IO_UNIVERSALRANK,IO_LOCALRANK,IO_JOBID,IO_DOMAINID) &
246  !$omp shared(IS,IE,JS,JE,EPS,UNDEF,Rdry,CPdry,PRC_myrank,IO_FID_LOG,IO_L,model_name,bulkflux, &
247  !$omp CPL_PHY_SFC_SKIN_itr_max,CPL_PHY_SFC_SKIN_dTS_max,CPL_PHY_SFC_SKIN_err_min,CPL_PHY_SFC_SKIN_res_min, &
248  !$omp calc_flag,dt,QVA,QVS,TMPA,TMPS,PRSA,RHOA,WA,UA,VA,LH,Z1,PBL, &
249  !$omp TG,PRSS,RHOS,TMPS1,WSTR,QVEF,Z0M,Z0H,Z0E,Rb,TC_dZ,ALBEDO,RFLXD, &
250  !$omp FracU10,FracT2,FracQ2, &
251  !$omp ZMFLX,XMFLX,YMFLX,SHFLX,LHFLX,QVFLX,GFLX,Ustar,Tstar,Qstar,Wstar,RLmo,U10,V10,T2,Q2) &
252 #else
253  !$omp default(shared) &
254 #endif
255  !$omp private(qdry,Rtot,flx_qv,redf,res,dts,olddts,emis,LWD,LWU,SWD,SWU,dres,oldres,dQVS, &
256  !$omp QVsat,dQVsat,dUstar,dTstar,dQstar,dWstar,dFracU10,dFracT2,dFracQ2, &
257  !$omp Uabs,dUabs,dRLmo,Ra,dRa,MFLUX)
258  do j = js, je
259  do i = is, ie
260  if ( calc_flag(i,j) ) then
261 
262 ! qdry = 1.0_RP - QVA(i,j)
263 ! Rtot = qdry * Rdry + QVA(i,j) * Rvap
264 
265  redf = 1.0_rp
266  oldres = huge(0.0_rp)
267  olddts = cpl_phy_sfc_skin_dts_max * dt
268 
269  ! modified Newton-Raphson method (Tomita 2009)
270  do n = 1, cpl_phy_sfc_skin_itr_max
271 
272  call qsat( tmps1(i,j), rhos(i,j), qvsat )
273  call qsat( tmps1(i,j)+dts0, rhos(i,j), dqvsat )
274 ! call qsat( TMPS1(i,j), PRSS(i,j), qdry, QVsat )
275 ! call qsat( TMPS1(i,j)+dTS0, PRSS(i,j), qdry, dQVsat )
276 
277  qvs(i,j) = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
278  + ( qvef(i,j) ) * qvsat
279  dqvs = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
280  + ( qvef(i,j) ) * dqvsat
281 
282  uabs = sqrt( wa(i,j)**2 + ua(i,j)**2 + va(i,j)**2 )
283 
284  call bulkflux( tmpa(i,j), tmps1(i,j), & ! [IN]
285  prsa(i,j), prss(i,j), & ! [IN]
286  qva(i,j), qvs(i,j), & ! [IN]
287  uabs, z1(i,j), pbl(i,j), & ! [IN]
288  z0m(i,j), z0h(i,j), z0e(i,j), & ! [IN]
289  ustar(i,j), tstar(i,j), qstar(i,j), & ! [OUT]
290  wstar(i,j), rlmo(i,j), ra, & ! [OUT]
291  fracu10(i,j), fract2(i,j), fracq2(i,j) ) ! [OUT]
292 
293  call bulkflux( tmpa(i,j), tmps1(i,j)+dts0, & ! [IN]
294  prsa(i,j), prss(i,j), & ! [IN]
295  qva(i,j), dqvs, & ! [IN]
296  uabs, z1(i,j), pbl(i,j), & ! [IN]
297  z0m(i,j), z0h(i,j), z0e(i,j), & ! [IN]
298  dustar, dtstar, dqstar, & ! [OUT]
299  dwstar, drlmo, dra, & ! [OUT] ! not used
300  dfracu10, dfract2, dfracq2 ) ! [OUT] ! not used
301 
302  emis = ( 1.0_rp - albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps1(i,j)**4
303 
304  lwd = rflxd(i,j,i_r_diffuse,i_r_ir)
305  lwu = rflxd(i,j,i_r_diffuse,i_r_ir) * albedo(i,j,i_r_diffuse,i_r_ir) + emis
306  swd = rflxd(i,j,i_r_direct ,i_r_nir) &
307  + rflxd(i,j,i_r_diffuse,i_r_nir) &
308  + rflxd(i,j,i_r_direct ,i_r_vis) &
309  + rflxd(i,j,i_r_diffuse,i_r_vis)
310  swu = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
311  + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
312  + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
313  + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)
314 
315  ! calculation for residual
316  flx_qv = min( - rhos(i,j) * ustar(i,j) * qstar(i,j) * ra / ( ra+rb(i,j) ), wstr(i,j)/real(dt,rp) )
317  res = swd - swu + lwd - lwu &
318  + cpdry * rhos(i,j) * ustar(i,j) * tstar(i,j) &
319  - lh(i,j) * flx_qv &
320  - tc_dz(i,j) * ( tmps1(i,j) - tg(i,j) )
321 
322  ! calculation for d(residual)/dTMPS
323  dres = -4.0_rp * emis / tmps1(i,j) &
324  + cpdry * rhos(i,j) * ( ustar(i,j)*(dtstar-tstar(i,j))/dts0 + tstar(i,j)*(dustar-ustar(i,j))/dts0 ) &
325  + lh(i,j) * rhos(i,j) * ( ustar(i,j)*(dqstar-qstar(i,j))/dts0 + qstar(i,j)*(dustar-ustar(i,j))/dts0 ) * ra / ( ra+rb(i,j) ) &
326  - tc_dz(i,j)
327 
328  ! convergence test with residual and error levels
329  if ( abs(res ) < cpl_phy_sfc_skin_res_min &
330  .OR. abs(res/dres) < cpl_phy_sfc_skin_err_min ) then
331  exit
332  endif
333 
334  ! calculate reduced factor
335  if ( dres < 0.0_rp ) then
336  if ( abs(res) > abs(oldres) ) then
337  redf = max( tfa*abs(redf), redf_min )
338  else
339  redf = min( tfb*abs(redf), redf_max )
340  endif
341  else
342  redf = -1.0_rp
343  endif
344 
345  ! estimate next surface temperature
346  dts = - redf * res / dres
347  dts = sign( min( abs(dts), abs(olddts) ), dts )
348  tmps1(i,j) = tmps1(i,j) + dts
349 
350  ! save residual in this step
351  oldres = res
352  olddts = dts
353  enddo
354 
355  ! update surface temperature with limitation
356  tmps1(i,j) = min( max( tmps1(i,j), &
357  tmps(i,j) - cpl_phy_sfc_skin_dts_max * real(dt,kind=rp) ), &
358  tmps(i,j) + cpl_phy_sfc_skin_dts_max * real(dt,kind=rp) )
359 
360  if ( n > cpl_phy_sfc_skin_itr_max ) then
361  ! surface temperature was not converged
362  log_warn("CPL_PHY_SFC_skin",*) 'surface tempearture was not converged. ', trim(model_name)
363  log_newline
364  log_info_cont('(A,I32)' ) 'number of i [no unit] :', i
365  log_info_cont('(A,I32)' ) 'number of j [no unit] :', j
366  log_newline
367  log_info_cont('(A,I32)' ) 'loop number [no unit] :', n
368  log_info_cont('(A,F32.16)') 'Residual [J/m2/s] :', res
369  log_info_cont('(A,F32.16)') 'delta Residual [J/m2/s] :', dres
370  log_newline
371  log_info_cont('(A,F32.16)') 'temperature [K] :', tmpa(i,j)
372  log_info_cont('(A,F32.16)') 'pressure [Pa] :', prsa(i,j)
373  log_info_cont('(A,F32.16)') 'velocity w [m/s] :', wa(i,j)
374  log_info_cont('(A,F32.16)') 'velocity u [m/s] :', ua(i,j)
375  log_info_cont('(A,F32.16)') 'velocity v [m/s] :', va(i,j)
376  log_info_cont('(A,F32.16)') 'absolute velocity [m/s] :', uabs
377  log_info_cont('(A,F32.16)') 'density [kg/m3] :', rhoa(i,j)
378  log_info_cont('(A,F32.16)') 'water vapor mass ratio [kg/kg] :', qva(i,j)
379  log_info_cont('(A,F32.16)') 'cell center height [m] :', z1(i,j)
380  log_info_cont('(A,F32.16)') 'atmospheric mixing layer height [m] :', pbl(i,j)
381  log_info_cont('(A,F32.16)') 'pressure at the surface [Pa] :', prss(i,j)
382  log_info_cont('(A,F32.16)') 'downward radiation (IR, direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_ir )
383  log_info_cont('(A,F32.16)') 'downward radiation (IR, diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_ir )
384  log_info_cont('(A,F32.16)') 'downward radiation (NIR,direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_nir)
385  log_info_cont('(A,F32.16)') 'downward radiation (NIR,diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_nir)
386  log_info_cont('(A,F32.16)') 'downward radiation (VIS,direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_vis)
387  log_info_cont('(A,F32.16)') 'downward radiation (VIS,diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_vis)
388  log_newline
389  log_info_cont('(A,F32.16)') 'soil temperature [K] :', tg(i,j)
390  log_info_cont('(A,F32.16)') 'soil water [kg/m2] :', wstr(i,j)
391  log_info_cont('(A,F32.16)') 'surface temperature [K] :', tmps(i,j)
392  log_info_cont('(A,F32.16)') 'surface density [kg/m3] :', rhos(i,j)
393  log_info_cont('(A,F32.16)') 'efficiency of evaporation [1] :', qvef(i,j)
394  log_info_cont('(A,F32.16)') 'surface albedo (IR, direct ) [1] :', albedo(i,j,i_r_direct ,i_r_ir )
395  log_info_cont('(A,F32.16)') 'surface albedo (IR, diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_ir )
396  log_info_cont('(A,F32.16)') 'surface albedo (NIR,direct ) [1] :', albedo(i,j,i_r_direct ,i_r_nir)
397  log_info_cont('(A,F32.16)') 'surface albedo (NIR,diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_nir)
398  log_info_cont('(A,F32.16)') 'surface albedo (VIS,direct ) [1] :', albedo(i,j,i_r_direct ,i_r_vis)
399  log_info_cont('(A,F32.16)') 'surface albedo (VIS,diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_vis)
400  log_info_cont('(A,F32.16)') 'latent heat [J/kg] :', lh(i,j)
401  log_info_cont('(A,F32.16)') 'stomata registance [1/s] :', rb(i,j)
402  log_info_cont('(A,F32.16)') 'thermal conductivity / depth [J/m2/s/K] :', tc_dz(i,j)
403  log_info_cont('(A,F32.16)') 'roughness length for momemtum [m] :', z0m(i,j)
404  log_info_cont('(A,F32.16)') 'roughness length for heat [m] :', z0h(i,j)
405  log_info_cont('(A,F32.16)') 'roughness length for vapor [m] :', z0e(i,j)
406  log_info_cont('(A,F32.16)') 'time step [s] :', dt
407  log_newline
408  log_info_cont('(A,F32.16)') 'friction velocity [m/s] :', ustar(i,j)
409  log_info_cont('(A,F32.16)') 'friction potential temperature [K] :', tstar(i,j)
410  log_info_cont('(A,F32.16)') 'friction water vapor mass ratio [kg/kg] :', qstar(i,j)
411  log_info_cont('(A,F32.16)') 'free convection velocity scale [m/s] :', wstar(i,j)
412  log_info_cont('(A,F32.16)') 'd(friction velocity) [m/s] :', dustar
413  log_info_cont('(A,F32.16)') 'd(friction potential temperature) [K] :', dtstar
414  log_info_cont('(A,F32.16)') 'd(friction water vapor mass ratio) [kg/kg] :', dqstar
415  log_info_cont('(A,F32.16)') 'd(free convection velocity scale) [m/s] :', dwstar
416  log_info_cont('(A,F32.16)') 'next surface temperature [K] :', tmps1(i,j)
417 
418  ! check NaN
419  if ( .NOT. ( res > -1.0_rp .OR. res < 1.0_rp ) ) then ! must be NaN
420  log_error("CPL_PHY_SFC_skin",*) 'NaN is detected for surface temperature. ', trim(model_name)
421  log_error_cont('(A,I32)' ) 'number of i [no unit] :', i
422  log_error_cont('(A,I32)' ) 'number of j [no unit] :', j
423  log_error_cont('(A,I32)' ) 'loop number [no unit] :', n
424  log_error_cont('(A,F32.16)') 'temperature [K] :', tmpa(i,j)
425  log_error_cont('(A,F32.16)') 'pressure [Pa] :', prsa(i,j)
426  log_error_cont('(A,F32.16)') 'velocity w [m/s] :', wa(i,j)
427  log_error_cont('(A,F32.16)') 'velocity u [m/s] :', ua(i,j)
428  log_error_cont('(A,F32.16)') 'velocity v [m/s] :', va(i,j)
429  log_error_cont('(A,F32.16)') 'absolute velocity [m/s] :', uabs
430  log_error_cont('(A,F32.16)') 'density [kg/m3] :', rhoa(i,j)
431  log_error_cont('(A,F32.16)') 'water vapor mass ratio [kg/kg] :', qva(i,j)
432  log_error_cont('(A,F32.16)') 'cell center height [m] :', z1(i,j)
433  log_error_cont('(A,F32.16)') 'atmospheric mixing layer height [m] :', pbl(i,j)
434  log_error_cont('(A,F32.16)') 'pressure at the surface [Pa] :', prss(i,j)
435  log_error_cont('(A,F32.16)') 'downward radiation (IR, direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_ir )
436  log_error_cont('(A,F32.16)') 'downward radiation (IR, diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_ir )
437  log_error_cont('(A,F32.16)') 'downward radiation (NIR,direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_nir)
438  log_error_cont('(A,F32.16)') 'downward radiation (NIR,diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_nir)
439  log_error_cont('(A,F32.16)') 'downward radiation (VIS,direct ) [J/m2/s] :', rflxd(i,j,i_r_direct ,i_r_vis)
440  log_error_cont('(A,F32.16)') 'downward radiation (VIS,diffuse) [J/m2/s] :', rflxd(i,j,i_r_diffuse,i_r_vis)
441  log_error_cont('(A,F32.16)') 'soil temperature [K] :', tg(i,j)
442  log_error_cont('(A,F32.16)') 'soil water [kg/m2] :', wstr(i,j)
443  log_error_cont('(A,F32.16)') 'surface temperature [K] :', tmps(i,j)
444  log_error_cont('(A,F32.16)') 'surface density [kg/m3] :', rhos(i,j)
445  log_error_cont('(A,F32.16)') 'efficiency of evaporation [1] :', qvef(i,j)
446  log_error_cont('(A,F32.16)') 'surface albedo (IR, direct ) [1] :', albedo(i,j,i_r_direct ,i_r_ir )
447  log_error_cont('(A,F32.16)') 'surface albedo (IR, diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_ir )
448  log_error_cont('(A,F32.16)') 'surface albedo (NIR,direct ) [1] :', albedo(i,j,i_r_direct ,i_r_nir)
449  log_error_cont('(A,F32.16)') 'surface albedo (NIR,diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_nir)
450  log_error_cont('(A,F32.16)') 'surface albedo (VIS,direct ) [1] :', albedo(i,j,i_r_direct ,i_r_vis)
451  log_error_cont('(A,F32.16)') 'surface albedo (VIS,diffuse) [1] :', albedo(i,j,i_r_diffuse,i_r_vis)
452  log_error_cont('(A,F32.16)') 'latent heat [J/kg] :', lh(i,j)
453  log_error_cont('(A,F32.16)') 'stomata registance [1/s] :', rb(i,j)
454  log_error_cont('(A,F32.16)') 'thermal conductivity / depth [J/m2/s/K] :', tc_dz(i,j)
455  log_error_cont('(A,F32.16)') 'roughness length for momemtum [m] :', z0m(i,j)
456  log_error_cont('(A,F32.16)') 'roughness length for heat [m] :', z0h(i,j)
457  log_error_cont('(A,F32.16)') 'roughness length for vapor [m] :', z0e(i,j)
458  log_error_cont('(A,F32.16)') 'time step [s] :', dt
459  log_error_cont('(A,F32.16)') 'friction velocity [m/s] :', ustar(i,j)
460  log_error_cont('(A,F32.16)') 'friction potential temperature [K] :', tstar(i,j)
461  log_error_cont('(A,F32.16)') 'friction water vapor mass ratio [kg/kg] :', qstar(i,j)
462  log_error_cont('(A,F32.16)') 'free convection velocity scale [m/s] :', wstar(i,j)
463  log_error_cont('(A,F32.16)') 'd(friction velocity) [m/s] :', dustar
464  log_error_cont('(A,F32.16)') 'd(friction potential temperature) [K] :', dtstar
465  log_error_cont('(A,F32.16)') 'd(friction water vapor mass ratio) [kg/kg] :', dqstar
466  log_error_cont('(A,F32.16)') 'd(free convection velocity scale) [m/s] :', dwstar
467  log_error_cont('(A,F32.16)') 'next surface temperature [K] :', tmps1(i,j)
468  call prc_abort
469  endif
470  endif
471 
472  ! calculate surface flux
473  tmps(i,j) = tmps1(i,j)
474 
475 ! qdry = 1.0_RP - QVA(i,j)
476  ! Rtot = qdry * Rdry + QVA(i,j) * Rvap
477 
478  call qsat( tmps(i,j), rhos(i,j), qvsat )
479 ! call qsat( TMPS(i,j), PRSS(i,j), qdry, QVsat )
480 
481  qvs(i,j) = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
482  + ( qvef(i,j) ) * qvsat
483 
484  call bulkflux( tmpa(i,j), tmps(i,j), & ! [IN]
485  prsa(i,j), prss(i,j), & ! [IN]
486  qva(i,j), qvs(i,j), & ! [IN]
487  uabs, z1(i,j), pbl(i,j), & ! [IN]
488  z0m(i,j), z0h(i,j), z0e(i,j), & ! [IN]
489  ustar(i,j), tstar(i,j), qstar(i,j), & ! [OUT]
490  wstar(i,j), rlmo(i,j), ra, & ! [OUT]
491  fracu10(i,j), fract2(i,j), fracq2(i,j) ) ! [OUT]
492 
493  if ( uabs < eps ) then
494  zmflx(i,j) = 0.0_rp
495  xmflx(i,j) = 0.0_rp
496  ymflx(i,j) = 0.0_rp
497  else
498  mflux = - rhos(i,j) * ustar(i,j)**2
499  zmflx(i,j) = mflux * wa(i,j) / uabs
500  xmflx(i,j) = mflux * ua(i,j) / uabs
501  ymflx(i,j) = mflux * va(i,j) / uabs
502  end if
503  shflx(i,j) = -rhos(i,j) * ustar(i,j) * tstar(i,j) * cpdry
504  qvflx(i,j) = min( - rhos(i,j) * ustar(i,j) * qstar(i,j) * ra / ( ra+rb(i,j) ), wstr(i,j)/real(dt,rp) )
505 
506  emis = ( 1.0_rp-albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps(i,j)**4
507 
508  lwd = rflxd(i,j,i_r_diffuse,i_r_ir)
509  lwu = rflxd(i,j,i_r_diffuse,i_r_ir) * albedo(i,j,i_r_diffuse,i_r_ir) + emis
510  swd = rflxd(i,j,i_r_direct ,i_r_nir) &
511  + rflxd(i,j,i_r_diffuse,i_r_nir) &
512  + rflxd(i,j,i_r_direct ,i_r_vis) &
513  + rflxd(i,j,i_r_diffuse,i_r_vis)
514  swu = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
515  + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
516  + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
517  + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)
518 
519  gflx(i,j) = tc_dz(i,j) * ( tmps(i,j) - tg(i,j) )
520 
521  lhflx(i,j) = qvflx(i,j) * lh(i,j)
522 
523 
524  ! calculation for residual
525  res = swd - swu + lwd - lwu - shflx(i,j) - lhflx(i,j) - gflx(i,j)
526 
527  ! put residual in ground heat flux
528  gflx(i,j) = gflx(i,j) + res
529 
530  else ! not calculate surface flux
531  zmflx(i,j) = undef
532  xmflx(i,j) = undef
533  ymflx(i,j) = undef
534  shflx(i,j) = undef
535  lhflx(i,j) = undef
536  qvflx(i,j) = undef
537  gflx(i,j) = undef
538  ustar(i,j) = undef
539  tstar(i,j) = undef
540  qstar(i,j) = undef
541  wstar(i,j) = undef
542  rlmo(i,j) = undef
543  u10(i,j) = undef
544  v10(i,j) = undef
545  t2(i,j) = undef
546  q2(i,j) = undef
547  endif
548  enddo
549  enddo
550 
551  call bulkflux_diagnose_surface( ia, is, ie, ja, js, je, &
552  ua(:,:), va(:,:), & ! (in)
553  tmpa(:,:), qva(:,:), & ! (in)
554  tmps(:,:), qvs(:,:), & ! (in)
555  z1(:,:), z0m(:,:), z0h(:,:), z0e(:,:), & ! (in)
556  u10(:,:), v10(:,:), t2(:,:), q2(:,:), & ! (out)
557  mask = calc_flag(:,:), & ! (in)
558  fracu10 = fracu10(:,:), & ! (in)
559  fract2 = fract2(:,:), & ! (in)
560  fracq2 = fracq2(:,:) ) ! (in)
561 
562  return
563  end subroutine cpl_phy_sfc_skin
564 
565 end module scale_cpl_phy_sfc_skin
scale_cpl_sfc_index::n_rad_dir
integer, parameter, public n_rad_dir
Definition: scale_cpl_sfc_index.F90:36
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_cpl_sfc_index::i_r_direct
integer, parameter, public i_r_direct
Definition: scale_cpl_sfc_index.F90:37
scale_cpl_sfc_index::i_r_diffuse
integer, parameter, public i_r_diffuse
Definition: scale_cpl_sfc_index.F90:38
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_cpl_sfc_index::i_r_ir
integer, parameter, public i_r_ir
Definition: scale_cpl_sfc_index.F90:29
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_bulkflux
module Surface bulk flux
Definition: scale_bulkflux.F90:12
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:90
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_io
module STDIO
Definition: scale_io.F90:10
scale_cpl_sfc_index::i_r_nir
integer, parameter, public i_r_nir
Definition: scale_cpl_sfc_index.F90:30
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_cpl_phy_sfc_skin
module coupler / physics / surface skin
Definition: scale_cpl_phy_sfc_skin.F90:12
scale_const::const_stb
real(rp), parameter, public const_stb
Stefan-Boltzman constant [W/m2/K4].
Definition: scale_const.F90:49
scale_cpl_phy_sfc_skin::cpl_phy_sfc_skin
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, WSTR, QVEF, ALBEDO, Rb, TC_dZ, Z0M, Z0H, Z0E, calc_flag, dt, model_name, TMPS, ZMFLX, XMFLX, YMFLX, SHFLX, LHFLX, QVFLX, GFLX, Ustar, Tstar, Qstar, Wstar, RLmo, U10, V10, T2, Q2)
Definition: scale_cpl_phy_sfc_skin.F90:114
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:90
scale_cpl_sfc_index
module coupler / surface-atmospehre
Definition: scale_cpl_sfc_index.F90:11
scale_cpl_sfc_index::i_r_vis
integer, parameter, public i_r_vis
Definition: scale_cpl_sfc_index.F90:31
scale_atmos_hydrometeor::lhf
real(rp), public lhf
latent heat of fusion for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:128
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_cpl_phy_sfc_skin::cpl_phy_sfc_skin_setup
subroutine, public cpl_phy_sfc_skin_setup
Setup.
Definition: scale_cpl_phy_sfc_skin.F90:56
scale_cpl_sfc_index::n_rad_rgn
integer, parameter, public n_rad_rgn
Definition: scale_cpl_sfc_index.F90:28
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_atmos_hydrometeor::cv_water
real(rp), public cv_water
CV for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:132
scale_atmos_hydrometeor::cv_ice
real(rp), public cv_ice
CV for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:134
scale_bulkflux::bulkflux
procedure(bc), pointer, public bulkflux
Definition: scale_bulkflux.F90:78
Generated by   doxygen 1.8.17