119 undef => const_undef, &
120 pre00 => const_pre00, &
121 tem00 => const_tem00, &
122 rdry => const_rdry, &
123 cpdry => const_cpdry, &
124 rvap => const_rvap, &
127 qsat => atmos_saturation_dens2qsat_all
130 atmos_hydrometeor_lhv, &
131 atmos_hydrometeor_lhs, &
137 bulkflux_diagnose_surface
140 integer,
intent(in) :: IA, IS, IE
141 integer,
intent(in) :: JA, JS, JE
142 real(RP),
intent(in) :: TMPA (IA,JA)
143 real(RP),
intent(in) :: PRSA (IA,JA)
144 real(RP),
intent(in) :: WA (IA,JA)
145 real(RP),
intent(in) :: UA (IA,JA)
146 real(RP),
intent(in) :: VA (IA,JA)
147 real(RP),
intent(in) :: RHOA (IA,JA)
148 real(RP),
intent(in) :: QVA (IA,JA)
149 real(RP),
intent(in) :: LH (IA,JA)
150 real(RP),
intent(in) :: Z1 (IA,JA)
151 real(RP),
intent(in) :: PBL (IA,JA)
152 real(RP),
intent(in) :: RHOS (IA,JA)
153 real(RP),
intent(in) :: PRSS (IA,JA)
154 real(RP),
intent(in) :: RFLXD (IA,JA,N_RAD_DIR,N_RAD_RGN)
155 real(RP),
intent(in) :: TG (IA,JA)
156 real(RP),
intent(in) :: WSTR (IA,JA)
157 real(RP),
intent(in) :: QVEF (IA,JA)
158 real(RP),
intent(in) :: ALBEDO (IA,JA,N_RAD_DIR,N_RAD_RGN)
159 real(RP),
intent(in) :: Rb (IA,JA)
160 real(RP),
intent(in) :: TC_dZ (IA,JA)
161 real(RP),
intent(in) :: Z0M (IA,JA)
162 real(RP),
intent(in) :: Z0H (IA,JA)
163 real(RP),
intent(in) :: Z0E (IA,JA)
164 logical,
intent(in) :: calc_flag(IA,JA)
165 real(DP),
intent(in) :: dt
166 character(len=*),
intent(in) :: model_name
168 real(RP),
intent(inout) :: TMPS (IA,JA)
170 real(RP),
intent(out) :: ZMFLX (IA,JA)
171 real(RP),
intent(out) :: XMFLX (IA,JA)
172 real(RP),
intent(out) :: YMFLX (IA,JA)
173 real(RP),
intent(out) :: SHFLX (IA,JA)
174 real(RP),
intent(out) :: LHFLX (IA,JA)
175 real(RP),
intent(out) :: QVFLX (IA,JA)
176 real(RP),
intent(out) :: GFLX (IA,JA)
177 real(RP),
intent(out) :: Ustar (IA,JA)
178 real(RP),
intent(out) :: Tstar (IA,JA)
179 real(RP),
intent(out) :: Qstar (IA,JA)
180 real(RP),
intent(out) :: Wstar (IA,JA)
181 real(RP),
intent(out) :: RLmo (IA,JA)
182 real(RP),
intent(out) :: U10 (IA,JA)
183 real(RP),
intent(out) :: V10 (IA,JA)
184 real(RP),
intent(out) :: T2 (IA,JA)
185 real(RP),
intent(out) :: Q2 (IA,JA)
187 real(RP),
parameter :: dTS0 = 1.0e-4_rp
188 real(RP),
parameter :: redf_min = 1.0e-2_rp
189 real(RP),
parameter :: redf_max = 1.0e+0_rp
190 real(RP),
parameter :: TFa = 0.5e+0_rp
191 real(RP),
parameter :: TFb = 1.1e+0_rp
193 real(RP) :: TMPS1(IA,JA)
214 real(RP) :: Uabs, dUabs
217 real(RP) :: QVsat, dQVsat
218 real(RP) :: QVS(IA,JA), dQVS
222 real(RP) :: FracU10(IA,JA), dFracU10
223 real(RP) :: FracT2 (IA,JA), dFracT2
224 real(RP) :: FracQ2 (IA,JA), dFracQ2
231 log_progress(*)
'coupler / physics / surface / SKIN'
237 tmps1(i,j) = tmps(i,j)
260 if ( calc_flag(i,j) )
then
266 oldres = huge(0.0_rp)
267 olddts = cpl_phy_sfc_skin_dts_max * dt
270 do n = 1, cpl_phy_sfc_skin_itr_max
272 call qsat( tmps1(i,j), rhos(i,j), qvsat )
273 call qsat( tmps1(i,j)+dts0, rhos(i,j), dqvsat )
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
282 uabs = sqrt( wa(i,j)**2 + ua(i,j)**2 + va(i,j)**2 )
284 call bulkflux( tmpa(i,j), tmps1(i,j), &
285 prsa(i,j), prss(i,j), &
286 qva(i,j), qvs(i,j), &
287 uabs, z1(i,j), pbl(i,j), &
288 z0m(i,j), z0h(i,j), z0e(i,j), &
289 ustar(i,j), tstar(i,j), qstar(i,j), &
290 wstar(i,j), rlmo(i,j), ra, &
291 fracu10(i,j), fract2(i,j), fracq2(i,j) )
293 call bulkflux( tmpa(i,j), tmps1(i,j)+dts0, &
294 prsa(i,j), prss(i,j), &
296 uabs, z1(i,j), pbl(i,j), &
297 z0m(i,j), z0h(i,j), z0e(i,j), &
298 dustar, dtstar, dqstar, &
299 dwstar, drlmo, dra, &
300 dfracu10, dfract2, dfracq2 )
302 emis = ( 1.0_rp - albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps1(i,j)**4
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)
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) &
320 - tc_dz(i,j) * ( tmps1(i,j) - tg(i,j) )
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) ) &
329 if ( abs(res ) < cpl_phy_sfc_skin_res_min &
330 .OR. abs(res/dres) < cpl_phy_sfc_skin_err_min )
then
335 if ( dres < 0.0_rp )
then
336 if ( abs(res) > abs(oldres) )
then
337 redf = max( tfa*abs(redf), redf_min )
339 redf = min( tfb*abs(redf), redf_max )
346 dts = - redf * res / dres
347 dts = sign( min( abs(dts), abs(olddts) ), dts )
348 tmps1(i,j) = tmps1(i,j) + dts
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) )
360 if ( n > cpl_phy_sfc_skin_itr_max )
then
362 log_warn(
"CPL_PHY_SFC_skin",*)
'surface tempearture was not converged. ', trim(model_name)
364 log_info_cont(
'(A,I32)' )
'number of i [no unit] :', i
365 log_info_cont(
'(A,I32)' )
'number of j [no unit] :', j
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
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)
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
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)
419 if ( .NOT. ( res > -1.0_rp .OR. res < 1.0_rp ) )
then
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)
473 tmps(i,j) = tmps1(i,j)
478 call qsat( tmps(i,j), rhos(i,j), qvsat )
481 qvs(i,j) = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
482 + ( qvef(i,j) ) * qvsat
484 call bulkflux( tmpa(i,j), tmps(i,j), &
485 prsa(i,j), prss(i,j), &
486 qva(i,j), qvs(i,j), &
487 uabs, z1(i,j), pbl(i,j), &
488 z0m(i,j), z0h(i,j), z0e(i,j), &
489 ustar(i,j), tstar(i,j), qstar(i,j), &
490 wstar(i,j), rlmo(i,j), ra, &
491 fracu10(i,j), fract2(i,j), fracq2(i,j) )
493 if ( uabs < eps )
then
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
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) )
506 emis = ( 1.0_rp-albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps(i,j)**4
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)
519 gflx(i,j) = tc_dz(i,j) * ( tmps(i,j) - tg(i,j) )
521 lhflx(i,j) = qvflx(i,j) * lh(i,j)
525 res = swd - swu + lwd - lwu - shflx(i,j) - lhflx(i,j) - gflx(i,j)
528 gflx(i,j) = gflx(i,j) + res
551 call bulkflux_diagnose_surface( ia, is, ie, ja, js, je, &
553 tmpa(:,:), qva(:,:), &
554 tmps(:,:), qvs(:,:), &
555 z1(:,:), z0m(:,:), z0h(:,:), z0e(:,:), &
556 u10(:,:), v10(:,:), t2(:,:), q2(:,:), &
557 mask = calc_flag(:,:), &
558 fracu10 = fracu10(:,:), &
559 fract2 = fract2(:,:), &
560 fracq2 = fracq2(:,:) )