source: issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90@ 27646

Last change on this file since 27646 was 27646, checked in by inwoo, 2 years ago

CHG: fix SMBsemic - consider QMR(heat from melting or refreezing) in semic..
File size: 10.8 KB
Line 
1subroutine run_semic_transient(nx, ntime, nloop, sf_in, rf_in, swd_in, lwd_in, wind_in, &
2 sp_in, rhoa_in, qq_in, tt_in, tsurf_in, qmr_in, &
3 tstic, &
4 hcrit, rcrit, &
5 mask, hice, hsnow, &
6 albedo, albedo_snow, &
7 alb_scheme, alb_smax, alb_smin, albi, albl, &
8 Tamp, &
9 tmin, tmax, tmid, mcrit, w_crit, tau_a, tau_f, afac, verbose, &
10 tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, alb_out, &
11 alb_snow_out,hsnow_out,hice_out,qmr_out) !{{{
12
13 use utils
14 use surface_physics
15 implicit none
16
17 ! declare surface physics class
18 type(surface_physics_class) :: surface
19 ! declare forcing class
20 !type(forc_class) :: forc
21 ! declare validation class
22 !type(vali_class) :: vali ! validation not needed here
23
24 integer, parameter:: dp=kind(0.d0) !< define precision (machine specific)
25 integer :: i, k, n
26 integer :: nnx, nny
27 integer :: year=0
28 integer :: day =1 !< not used value.
29 integer, intent(in) :: nx, ntime, nloop ! number of grid / number of times
30 logical, intent(in) :: verbose ! verbosity
31 logical :: debug=.false.
32
33 ! forcing data
34 ! input argument format array size (nx, ntime)...
35 double precision, intent(in), dimension(nx):: sf_in ! snow fall rate [m/s]
36 double precision, intent(in), dimension(nx):: rf_in ! rain fall rate [m/s]
37 double precision, intent(in), dimension(nx):: swd_in ! downwelling shortwave radiation [W/m2]
38 double precision, intent(in), dimension(nx):: lwd_in ! downwelling longwave radiation [W/m2]
39 double precision, intent(in), dimension(nx):: wind_in ! surface wind speed [m/s]
40 double precision, intent(in), dimension(nx):: sp_in ! surface pressure [Pa]
41 double precision, intent(in), dimension(nx):: rhoa_in ! air density [kg/m3]
42 double precision, intent(in), dimension(nx):: qq_in ! air specific humidity [kg/kg]
43 double precision, intent(in), dimension(nx):: tt_in ! air temperature [K]
44
45 ! input data
46 double precision, intent(in) :: tstic ! time step from ISSM [sec].
47
48 ! output data
49 ! Ice surface Temperature [K]
50 double precision, intent(out), dimension(nx):: tsurf_out
51 ! surface mass balance=(Accu-Melt) [m/s]
52 double precision, intent(out), dimension(nx):: smb_out
53 double precision, intent(out), dimension(nx):: smbi_out ! SMB ice [water equivalent m/s]
54 double precision, intent(out), dimension(nx):: smbs_out ! SMB snow [water equivalent m/s]
55 double precision, intent(out), dimension(nx):: saccu_out ! accumulation [m/s]
56 double precision, intent(out), dimension(nx):: smelt_out ! ablation [m/s]
57 double precision, intent(out), dimension(nx):: alb_out ! grid-averaged albedo [no unit]
58 double precision, intent(out), dimension(nx):: alb_snow_out
59 double precision, intent(out), dimension(nx):: hice_out
60 double precision, intent(out), dimension(nx):: hsnow_out
61 double precision, intent(out), dimension(nx):: qmr_out
62
63 double precision :: total_time, start, finish
64
65 ! set parameters
66 !character (len=256) :: name ! not used(?)
67 !character (len=256) :: boundary(30) ! not used(?)
68 !character (len=256), intent(in) :: alb_scheme !< name of albedo scheme
69 integer, intent(in) :: alb_scheme
70 !integer :: n_ksub
71 !double precision :: ceff !< surface heat heat capacity of snow/ice [J/W m2]
72 double precision, intent(in), dimension(nx):: albedo
73 double precision, intent(in), dimension(nx):: albedo_snow !< spatial..
74 double precision, intent(in), dimension(nx):: hsnow
75 double precision, intent(in), dimension(nx):: hice
76 double precision, intent(in), dimension(nx):: tsurf_in !< input temperature [K]
77 double precision, intent(in), dimension(nx):: qmr_in
78 double precision, intent(in), dimension(nx):: mask
79
80 double precision, intent(in) :: albi
81 double precision, intent(in) :: albl
82 double precision, intent(in) :: alb_smax
83 double precision, intent(in) :: alb_smin
84 double precision, intent(in) :: hcrit !< critical snow height for which grid cell is 50% snow covered [m]
85 double precision, intent(in) :: rcrit !< critical snow height fro which refreezing fraction is 50% [m]
86 double precision, intent(in) :: Tamp
87 !double precision :: csh
88 !double precision :: clh
89 double precision, intent(in) :: tmin
90 double precision, intent(in) :: tmax
91 !double precision :: tsticsub
92 ! parameters for isba albedo scheme.
93 double precision, intent(in) :: tau_a !< critical liquide water concent for "isba" albedo scheme [kg/m2]
94 double precision, intent(in) :: tau_f
95 double precision, intent(in) :: w_crit
96 double precision, intent(in) :: mcrit
97 double precision, intent(in) :: afac !< param
98 double precision, intent(in) :: tmid !< param for "alex" albedo parameterization [K]
99
100 if (debug) then
101 print*,' ntime: ', ntime
102 print*,' nx : ', nx
103 end if
104
105 ! set vector length
106 surface%par%nx = nx
107
108 ! FIXME should be user input
109 !boundary = "" "" ""
110 if (debug) then
111 print*, "run_semic_transient: initialize parameters."
112 end if
113 surface%par%tstic = tstic !< time step [s]
114 surface%par%ceff= 2.0e6_dp !< surface heat capacity of snow/ice [J/K/m2]
115 surface%par%csh = 2.0e-3_dp !< turbulent heat exchange coefficient
116 surface%par%clh = 5.0e-4_dp !< turbulent heat exchange coefficient [no unit]
117 surface%par%alb_smax = alb_smax !0.79_dp !< max snow albedo
118 surface%par%alb_smin = alb_smin !0.6_dp !< min snow albedo
119 surface%par%albi = albi ! 0.41_dp !< albedo for ice
120 surface%par%albl = albl ! 0.07_dp !< albedo for land
121 surface%par%tmin = tmin ! -999_dp
122 surface%par%tmax = tmax ! 273.15_dp
123 surface%par%hcrit = hcrit !0.028_dp !< critical snow height for which grid cell is 50 % snow colvered [m]
124 surface%par%rcrit = rcrit !0.85_dp !< refreezing fraction is 50% [m]
125 surface%par%amp = Tamp !3.0_dp !< amplitude of diurnal cycle [K]
126 if (alb_scheme == 0) then
127 surface%par%alb_scheme="none"
128 else if (alb_scheme == 1) then
129 surface%par%alb_scheme = "slater"
130 else if (alb_scheme == 2) then
131 surface%par%alb_scheme = "denby"
132 else if (alb_scheme == 3) then
133 surface%par%alb_scheme = "isba"
134 else
135 print*, "ERROR: current albedo scheme is not available."
136 call exit(1)
137 end if
138 surface%par%tau_a = tau_a !0.008_dp
139 surface%par%tau_f = tau_f !0.24_dp
140 surface%par%w_crit = w_crit !15.0_dp ! default value
141 surface%par%mcrit = mcrit !6.0e-8_dp
142 surface%par%n_ksub = 3 ! sub ...
143 ! snow albedo of alex
144 surface%par%afac = afac
145 surface%par%tmid = tmid
146
147 ! initialize sub-daily time step tsticsub
148 surface%par%tsticsub = surface%par%tstic / dble(surface%par%n_ksub)
149
150 ! allocate necessary arrays for surface_physics module
151 call surface_alloc(surface%now,surface%par%nx)
152
153 ! initialise prognostic variables
154 if (debug) then
155 print*,"run_semic_transient: initialize variables."
156 end if
157 ! these values will be updated through "surface_energy_and_mass_balance" function.
158 surface%now%mask (:) = mask (:) ! 2.0_dp !loi_mask(:nx)
159 if (debug) then
160 print*,"run_semic_transient: initialize variables: mask"
161 end if
162 surface%now%hsnow (:) = hsnow (:) ! initial snow height...
163 surface%now%hice (:) = hice (:) ! initial ice height..
164 surface%now%tsurf (:) = tsurf_in (:) !< initial ice surface temperature
165 surface%now%alb (:) = albedo (:) !< initial albedo for energy balance.
166 surface%now%alb_snow(:) = albedo_snow(:) !< initial albedo for ISBA albedo method.
167 if (debug) then
168 print*,"run_semic_transient: initialize variables. DONE."
169 end if
170
171 if (debug) then
172 !print*, "====== global variable =========="
173 !print*, "nloop :", nloop
174 !print*, "nx :", surface%par%nx
175 !print*, "====== parameters ======"
176 !print*, "csh :", surface%par%csh
177 !print*, "clh :", surface%par%clh
178 !print*, "albeo scheme :", surface%par%alb_scheme
179 !print*, "albeo ice :", surface%par%albi
180 !print*, "tstic :", surface%par%tstic
181 !print*, "tsticsub :", surface%par%tsticsub
182 !print*, "n_ksub :", surface%par%n_ksub
183 print*, "====== inputs ========="
184 print*, "hsnow :", hsnow
185 print*, "====== state variables ======"
186 print*, "hsnow :", surface%now%hsnow
187 print*, "hice :", surface%now%hice
188 print*, "albeo :", surface%now%alb
189 print*, "albeo snow :", surface%now%alb_snow
190 print*, "mask :", surface%now%mask
191 print*, "tsurf :", surface%now%tsurf
192 print*, "sf :", sf_in
193 end if
194
195 ! define boundary conditions (not used, here!)
196 call surface_boundary_define(surface%bnd,surface%par%boundary)
197 !call print_boundary_opt(surface%bnd)
198
199 ! input with single value
200 do k =1,nloop
201 do i =1,ntime
202 if (debug) then
203 print*,"run_semic_transient: forcing data: ntime = ", i
204 end if
205 surface%now%sf = sf_in !(:,i)
206 surface%now%rf = rf_in !(:,i)
207 surface%now%sp = sp_in !(:,i)
208 surface%now%lwd = lwd_in !(:,i)
209 surface%now%swd = swd_in !(:,i)
210 surface%now%wind = wind_in!(:,i)
211 surface%now%rhoa = rhoa_in!(:,i)
212 surface%now%t2m = tt_in !(:,i)
213 surface%now%qq = qq_in !(:,i)
214 ! qmr_res is used to "energy_balance" in semic.
215 surface%now%qmr_res = qmr_in
216
217 ! calculate prognostic and diagnsotic variables
218 call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,i,year)
219
220 if (debug) then
221 print*,"done..."
222 end if
223 if (k == nloop) then
224 tsurf_out = surface%now%tsurf
225 ! melt - potential surface melt [m/s]
226 ! smb = SMB_ice + SMB_snow
227 ! smbi - SMB_ice (water equivalent m/sec)
228 ! smbs - SMB_snow (water equivalent m/sec)
229 smb_out =surface%now%smb ! smb = smb_snow + smb_ice
230 smbi_out =surface%now%smb_ice ! Csi (snow>ice) - melted_ice + refrezon_rain.
231 smbs_out =surface%now%smb_snow ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow
232 saccu_out =surface%now%acc ! acc = snowfall - sublimiation - refreezing
233 smelt_out =surface%now%melt ! potential surface melt = melt_ice + melt_snow
234 alb_out =surface%now%alb
235 alb_snow_out =surface%now%alb_snow
236 hsnow_out =surface%now%hsnow
237 hice_out =surface%now%hice
238 qmr_out =surface%now%qmr_res
239 end if
240 end do
241 end do
242
243 ! de-allocate surface_physics arrays
244 call surface_dealloc(surface%now)
245
246end subroutine run_semic_transient ! }}}
Note: See TracBrowser for help on using the repository browser.