source: LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ice_optics_fu.F90 @ 5408

Last change on this file since 5408 was 4773, checked in by idelkadi, 12 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 4.8 KB
Line 
1! radiation_ice_optics_fu.F90 - Fu's scheme for ice optical properties
2!
3! (C) Copyright 2014- ECMWF.
4!
5! This software is licensed under the terms of the Apache Licence Version 2.0
6! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7!
8! In applying this licence, ECMWF does not waive the privileges and immunities
9! granted to it by virtue of its status as an intergovernmental organisation
10! nor does it submit to any jurisdiction.
11!
12! Author:  Robin Hogan
13! Email:   r.j.hogan@ecmwf.int
14!
15! Modifications
16!   2020-08-10  R. Hogan  Bounded re to be <= 100um and g to be < 1.0
17
18module radiation_ice_optics_fu
19
20  use parkind1, only : jprb
21
22  implicit none
23  public
24
25  ! The number of ice coefficients depends on the parameterization
26  integer, parameter :: NIceOpticsCoeffsFuSW  = 10
27  integer, parameter :: NIceOpticsCoeffsFuLW  = 11
28
29  ! Limits based on the range of validity of the parameterizations
30  real(jprb), parameter :: MaxAsymmetryFactor = 1.0_jprb - 10.0_jprb*epsilon(1.0_jprb)
31  real(jprb), parameter :: MaxEffectiveRadius = 100.0e-6_jprb ! metres
32
33contains
34
35  !---------------------------------------------------------------------
36  ! Compute shortwave ice-particle scattering properties using Fu
37  ! (1996) parameterization.  The asymmetry factor in band 14 goes
38  ! larger than one for re > 100.8 um, so we cap re at 100 um.
39  ! Asymmetry factor is capped at just less than 1 because if it is
40  ! exactly 1 then delta-Eddington scaling leads to a zero scattering
41  ! optical depth and then division by zero.
42  subroutine calc_ice_optics_fu_sw(nb, coeff, ice_wp, &
43       &  re, od, scat_od, g)
44
45    !use yomhook,  only : lhook, dr_hook, jphook
46
47    ! Number of bands
48    integer, intent(in)  :: nb
49    ! Coefficients read from a data file
50    real(jprb), intent(in) :: coeff(:,:)
51    ! Ice water path (kg m-2)
52    real(jprb), intent(in) :: ice_wp
53    ! Effective radius (m)
54    real(jprb), intent(in) :: re
55    ! Total optical depth, scattering optical depth and asymmetry factor
56    real(jprb), intent(out) :: od(nb), scat_od(nb), g(nb)
57
58    ! Fu's effective diameter (microns) and its inverse
59    real(jprb) :: de_um, inv_de_um
60    ! Ice water path in g m-2
61    real (jprb) :: iwp_gm_2
62
63    integer :: jb
64    !real(jphook) :: hook_handle
65
66    !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_sw',0,hook_handle)
67
68    ! Convert to effective diameter using the relationship in the IFS
69    de_um     = min(re, MaxEffectiveRadius) * (1.0e6_jprb / 0.64952_jprb)
70    inv_de_um = 1.0_jprb / de_um
71    iwp_gm_2  = ice_wp * 1000.0_jprb
72
73! Added for DWD (2020)
74!NEC$ shortloop
75    do jb = 1, nb
76      od(jb) = iwp_gm_2 * (coeff(jb,1) + coeff(jb,2) * inv_de_um)
77      scat_od(jb) = od(jb) * (1.0_jprb - (coeff(jb,3) + de_um*(coeff(jb,4) &
78         &  + de_um*(coeff(jb,5) + de_um*coeff(jb,6)))))
79      g(jb) = min(coeff(jb,7) + de_um*(coeff(jb,8) &
80         &  + de_um*(coeff(jb,9) + de_um*coeff(jb,10))), &
81         &  MaxAsymmetryFactor)
82    end do
83
84    !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_sw',1,hook_handle)
85
86  end subroutine calc_ice_optics_fu_sw
87
88
89  !---------------------------------------------------------------------
90  ! Compute longwave ice-particle scattering properties using Fu et
91  ! al. (1998) parameterization
92  subroutine calc_ice_optics_fu_lw(nb, coeff, ice_wp, &
93       &  re, od, scat_od, g)
94
95    !use yomhook,  only : lhook, dr_hook, jphook
96
97    ! Number of bands
98    integer, intent(in)  :: nb
99    ! Coefficients read from a data file
100    real(jprb), intent(in) :: coeff(:,:)
101    ! Ice water path (kg m-2)
102    real(jprb), intent(in) :: ice_wp
103    ! Effective radius (m)
104    real(jprb), intent(in) :: re
105    ! Total optical depth, scattering optical depth and asymmetry factor
106    real(jprb), intent(out) :: od(nb), scat_od(nb), g(nb)
107
108    ! Fu's effective diameter (microns) and its inverse
109    real(jprb) :: de_um, inv_de_um
110    ! Ice water path in g m-2
111    real (jprb) :: iwp_gm_2
112
113    integer :: jb
114    !real(jphook) :: hook_handle
115
116    !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_lw',0,hook_handle)
117
118    ! Convert to effective diameter using the relationship in the IFS
119    de_um = min(re, MaxEffectiveRadius) * (1.0e6_jprb / 0.64952_jprb)
120
121    inv_de_um = 1.0_jprb / de_um
122    iwp_gm_2  = ice_wp * 1000.0_jprb
123
124! Added for DWD (2020)
125!NEC$ shortloop
126    do jb = 1, nb
127      od(jb) = iwp_gm_2 * (coeff(jb,1) + inv_de_um*(coeff(jb,2) &
128         &  + inv_de_um*coeff(jb,3)))
129      scat_od(jb) = od(jb) - iwp_gm_2*inv_de_um*(coeff(jb,4) + de_um*(coeff(jb,5) &
130         &  + de_um*(coeff(jb,6) + de_um*coeff(jb,7))))
131      g(jb) = min(coeff(jb,8) + de_um*(coeff(jb,9) &
132         &  + de_um*(coeff(jb,10) + de_um*coeff(jb,11))), &
133         &  MaxAsymmetryFactor)
134    end do
135
136    !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_lw',1,hook_handle)
137
138  end subroutine calc_ice_optics_fu_lw
139
140end module radiation_ice_optics_fu
Note: See TracBrowser for help on using the repository browser.