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