84 module module_radiation_surface
90 use module_iounitdef
, only : niradsf
97 character(40),
parameter :: &
98 & VTAGSFC=
'NCEP-Radiation_surface v5.1 Nov 2012 '
106 integer,
parameter,
public ::
imxems = 360
109 integer,
parameter,
public ::
jmxems = 180
111 real (kind=kind_phys),
parameter :: f_zero = 0.0
112 real (kind=kind_phys),
parameter :: f_one = 1.0
113 real (kind=kind_phys),
parameter :: rad2dg= 180.0 /
con_pi
168 integer,
intent(in) :: me
175 logical :: file_exist
176 character :: cline*80
180 if ( me == 0 ) print *, vtagsfc
190 print *,
' - Using climatology surface albedo scheme for sw'
196 print *,
' - Using MODIS based land surface albedo for sw'
200 print *,
' !! ERROR in Albedo Scheme Setting, IALB=',
ialbflg
213 print *,
' - Using Fixed Surface Emissivity = 1.0 for lw'
216 elseif (
iemslw == 1 )
then
219 if ( .not.
allocated(
idxems) )
then
227 if ( .not. file_exist )
then
229 print *,
' - Using Varying Surface Emissivity for lw'
230 print *,
' Requested data file "',
semis_file,
'" not found!'
231 print *,
' Change to fixed surface emissivity = 1.0 !'
237 open (niradsf,file=
semis_file,form=
'formatted',status=
'old')
240 read (niradsf,12) cline
247 print *,
' - Using Varying Surface Emissivity for lw'
260 print *,
' !! ERROR in Emissivity Scheme Setting, IEMS=',
iemsflg
309 & ( slmsk,snowf,sncovr,snoalb,zorlf,coszf,tsknf,tairf,hprif, &
310 & alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc, &
375 integer,
intent(in) :: IMAX
377 real (kind=kind_phys),
dimension(:),
intent(in) :: &
378 & slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif, &
379 & alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, &
383 real (kind=kind_phys),
dimension(IMAX,NF_ALBD),
intent(out) :: &
388 real (kind=kind_phys) :: asnvb, asnnb, asnvd, asnnd, asevb &
389 &, asenb, asevd, asend, fsno, fsea, rfcs, rfcw, flnd &
390 &, asnow, argh, hrgh, fsno0, fsno1, flnd0, fsea0, csnow &
391 &, a1, a2, b1, b2, b3, ab1bm, ab2bm
393 real (kind=kind_phys) ffw, dtgd
409 asnow = 0.02*snowf(i)
410 argh = min(0.50, max(.025, 0.01*zorlf(i)))
411 hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) )
412 fsno0 = asnow / (argh + asnow) * hrgh
413 if (nint(slmsk(i))==0 .and. tsknf(i)>
con_tice) fsno0 = f_zero
414 fsno1 = f_one - fsno0
415 flnd0 = min(f_one, facsf(i)+facwf(i))
416 fsea0 = max(f_zero, f_one-flnd0)
423 if (tsknf(i) >= 271.5)
then
426 elseif (tsknf(i) < 271.1)
then
430 a1 = (tsknf(i) - 271.1)**2
432 asend = 0.65 - 3.6875*a1
437 if (nint(slmsk(i)) == 2)
then
438 ffw = f_one - fice(i)
439 if (ffw < f_one)
then
440 dtgd = max(f_zero, min(5.0, (
con_ttp-tisfc(i)) ))
447 asnvd = (0.70 + b1) * fice(i) + b3
448 asnnd = (0.60 + b1) * fice(i) + b3
449 asevd = 0.70 * fice(i) + b3
450 asend = 0.60 * fice(i) + b3
458 if (coszf(i) < 0.5)
then
459 csnow = 0.5 * (3.0 / (f_one+4.0*coszf(i)) - f_one)
460 asnvb = min( 0.98, asnvd+(1.0-asnvd)*csnow )
461 asnnb = min( 0.98, asnnd+(1.0-asnnd)*csnow )
469 if (coszf(i) > 0.0001)
then
472 rfcs = 2.14 / (f_one + 1.48*coszf(i))
476 asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) &
477 & + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
478 & * (coszf(i)-f_one))
491 a1 = alvsf(i) * facsf(i)
492 b1 = alvwf(i) * facwf(i)
493 a2 = alnsf(i) * facsf(i)
494 b2 = alnwf(i) * facwf(i)
495 ab1bm = a1*rfcs + b1*rfcw
496 ab2bm = a2*rfcs + b2*rfcw
497 sfcalb(i,1) = min(0.99, ab2bm) *flnd + asenb*fsea + asnnb*fsno
498 sfcalb(i,2) = (a2 + b2) * 0.96 *flnd + asend*fsea + asnnd*fsno
499 sfcalb(i,3) = min(0.99, ab1bm) *flnd + asevb*fsea + asnvb*fsno
500 sfcalb(i,4) = (a1 + b1) * 0.96 *flnd + asevd*fsea + asnvd*fsno
514 if (nint(slmsk(i))==0 .and. tsknf(i)>
con_tice) fsno0 = f_zero
516 if (nint(slmsk(i)) == 2)
then
517 asnow = 0.02*snowf(i)
518 argh = min(0.50, max(.025, 0.01*zorlf(i)))
519 hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) )
520 fsno0 = asnow / (argh + asnow) * hrgh
523 fsno1 = f_one - fsno0
524 flnd0 = min(f_one, facsf(i)+facwf(i))
525 fsea0 = max(f_zero, f_one-flnd0)
532 if (tsknf(i) >= 271.5)
then
535 elseif (tsknf(i) < 271.1)
then
539 a1 = (tsknf(i) - 271.1)**2
541 asend = 0.65 - 3.6875*a1
547 if (nint(slmsk(i)) == 2)
then
548 ffw = f_one - fice(i)
549 if (ffw < f_one)
then
550 dtgd = max(f_zero, min(5.0, (
con_ttp-tisfc(i)) ))
557 asnvd = (0.70 + b1) * fice(i) + b3
558 asnnd = (0.60 + b1) * fice(i) + b3
559 asevd = 0.70 * fice(i) + b3
560 asend = 0.60 * fice(i) + b3
568 if (nint(slmsk(i)) == 2)
then
569 if (coszf(i) < 0.5)
then
570 csnow = 0.5 * (3.0 / (f_one+4.0*coszf(i)) - f_one)
571 asnvb = min( 0.98, asnvd+(f_one-asnvd)*csnow )
572 asnnb = min( 0.98, asnnd+(f_one-asnnd)*csnow )
585 if (coszf(i) > 0.0001)
then
589 rfcs = 1.775/(1.0+1.55*coszf(i))
592 asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) &
593 & + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
594 & * (coszf(i)-f_one))
606 ab1bm = min(0.99, alnsf(i)*rfcs)
607 ab2bm = min(0.99, alvsf(i)*rfcs)
608 sfcalb(i,1) = ab1bm *flnd + asenb*fsea + asnnb*fsno
609 sfcalb(i,2) = alnwf(i) *flnd + asend*fsea + asnnd*fsno
610 sfcalb(i,3) = ab2bm *flnd + asevb*fsea + asnvb*fsno
611 sfcalb(i,4) = alvwf(i) *flnd + asevd*fsea + asnvd*fsno
641 & ( xlon,xlat,slmsk,snowf,sncovr,zorlf,tsknf,tairf,hprif, &
688 integer,
intent(in) :: IMAX
690 real (kind=kind_phys),
dimension(:),
intent(in) :: &
691 & xlon,xlat, slmsk, snowf,sncovr, zorlf, tsknf, tairf, hprif
694 real (kind=kind_phys),
dimension(:),
intent(out) :: sfcemis
697 integer :: i, i1, i2, j1, j2, idx
699 real (kind=kind_phys) :: dltg, hdlt, tmp1, tmp2, &
700 & asnow, argh, hrgh, fsno, fsno0, fsno1
706 real (kind=kind_phys) :: emsref(8)
707 data emsref / 0.97, 0.95, 0.94, 0.90, 0.93, 0.96, 0.96, 0.99 /
720 dltg = 360.0 / float(
imxems)
727 lab_do_imax :
do i = 1, imax
729 if ( nint(slmsk(i)) == 0 )
then
731 sfcemis(i) = emsref(1)
733 else if ( nint(slmsk(i)) == 2 )
then
735 sfcemis(i) = emsref(7)
743 tmp1 = xlon(i) * rad2dg
744 if (tmp1 < f_zero) tmp1 = tmp1 + 360.0
746 lab_do_imxems :
do i1 = 1,
imxems
747 tmp2 = dltg * (i1 - 1) + hdlt
749 if (abs(tmp1-tmp2) <= hdlt)
then
756 tmp1 = xlat(i) * rad2dg
759 lab_do_jmxems :
do j1 = 1,
jmxems
760 tmp2 = 90.0 - dltg * (j1 - 1)
762 if (abs(tmp1-tmp2) <= hdlt)
then
769 idx = max( 2,
idxems(i2,j2) )
770 if ( idx >= 7 ) idx = 2
771 sfcemis(i) = emsref(idx)
777 if (
ialbflg==1 .and. nint(slmsk(i))==1 )
then
780 fsno1 = f_one - fsno0
781 sfcemis(i) = sfcemis(i)*fsno1 + emsref(8)*fsno0
784 if ( snowf(i) > f_zero )
then
785 asnow = 0.02*snowf(i)
786 argh = min(0.50, max(.025, 0.01*zorlf(i)))
787 hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) )
788 fsno0 = asnow / (argh + asnow) * hrgh
789 if (nint(slmsk(i)) == 0 .and. tsknf(i) > 271.2) &
791 fsno1 = f_one - fsno0
792 sfcemis(i) = sfcemis(i)*fsno1 + emsref(8)*fsno0
812 end module module_radiation_surface
real(kind=kind_phys), parameter con_pi
pi
real(kind=kind_phys), parameter con_tice
temp freezing sea (K)
subroutine, public setemis(xlon, xlat, slmsk, snowf, sncovr, zorlf, tsknf, tairf, hprif, IMAX, sfcemis )
This subroutine computes surface emissivity for LW radiation.
integer, save ialbflg
surface albedo scheme control flag =0:vegetation type based climatological albedo scheme =1:seaso...
integer, dimension(:,:), allocatable idxems
global surface emissivity index array
subroutine, public sfc_init(me)
This subroutine is the initialization program for surface radiation related quantities (albedo...
character, save semis_file
external sfc emissivity data table: sfc_emissivity_idx.txt
integer, parameter, public nf_albd
num of sfc albedo components
real(kind=kind_phys), parameter con_t0c
temp at 0C (K)
integer, save iemsflg
surface emissivity scheme control flag =0:black-body surface emissivity(=1.0) =1:vegetation type ...
subroutine, public setalb(slmsk, snowf, sncovr, snoalb, zorlf, coszf, tsknf, tairf, hprif, alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, IMAX, sfcalb )
This subroutine computes four components of surface albedos (i.e., vis-nir, direct-diffused) accordin...
integer iemslw
global surface emissivity contrl flag set up in 'sfc_init'
integer, parameter, public imxems
num of longitude points in global emis-type map
real(kind=kind_phys), parameter con_ttp
temp at H2O 3pt (K)
integer, parameter, public jmxems
num of latitude points in global emis-type map