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

Last change on this file since 27773 was 27773, checked in by inwoo, 21 months ago

CHG: Update SEMIC to enable the output request of SmbRefreeze value. Also, synchronize the SMBsemic class for MATLAB and PYTHON.
File size: 11.0 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, refr_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):: refr_out ! freezing [m/s]
58 double precision, intent(out), dimension(nx):: alb_out ! grid-averaged albedo [no unit]
59 double precision, intent(out), dimension(nx):: alb_snow_out
60 double precision, intent(out), dimension(nx):: hice_out
61 double precision, intent(out), dimension(nx):: hsnow_out
62 double precision, intent(out), dimension(nx):: qmr_out
63
64 double precision :: total_time, start, finish
65
66 ! set parameters
67 !character (len=256) :: name ! not used(?)
68 !character (len=256) :: boundary(30) ! not used(?)
69 !character (len=256), intent(in) :: alb_scheme !< name of albedo scheme
70 integer, intent(in) :: alb_scheme
71 !integer :: n_ksub
72 !double precision :: ceff !< surface heat heat capacity of snow/ice [J/W m2]
73 double precision, intent(in), dimension(nx):: albedo
74 double precision, intent(in), dimension(nx):: albedo_snow !< spatial..
75 double precision, intent(in), dimension(nx):: hsnow
76 double precision, intent(in), dimension(nx):: hice
77 double precision, intent(in), dimension(nx):: tsurf_in !< input temperature [K]
78 double precision, intent(in), dimension(nx):: qmr_in
79 double precision, intent(in), dimension(nx):: mask
80
81 double precision, intent(in) :: albi
82 double precision, intent(in) :: albl
83 double precision, intent(in) :: alb_smax
84 double precision, intent(in) :: alb_smin
85 double precision, intent(in) :: hcrit !< critical snow height for which grid cell is 50% snow covered [m]
86 double precision, intent(in) :: rcrit !< critical snow height fro which refreezing fraction is 50% [m]
87 double precision, intent(in) :: Tamp
88 !double precision :: csh
89 !double precision :: clh
90 double precision, intent(in) :: tmin
91 double precision, intent(in) :: tmax
92 !double precision :: tsticsub
93 ! parameters for isba albedo scheme.
94 double precision, intent(in) :: tau_a !< critical liquide water concent for "isba" albedo scheme [kg/m2]
95 double precision, intent(in) :: tau_f
96 double precision, intent(in) :: w_crit
97 double precision, intent(in) :: mcrit
98 double precision, intent(in) :: afac !< param
99 double precision, intent(in) :: tmid !< param for "alex" albedo parameterization [K]
100
101 if (debug) then
102 print*,' ntime: ', ntime
103 print*,' nx : ', nx
104 end if
105
106 ! set vector length
107 surface%par%nx = nx
108
109 ! FIXME should be user input
110 !boundary = "" "" ""
111 if (debug) then
112 print*, "run_semic_transient: initialize parameters."
113 end if
114 surface%par%tstic = tstic !< time step [s]
115 surface%par%ceff= 2.0e6_dp !< surface heat capacity of snow/ice [J/K/m2]
116 surface%par%csh = 2.0e-3_dp !< turbulent heat exchange coefficient
117 surface%par%clh = 5.0e-4_dp !< turbulent heat exchange coefficient [no unit]
118 surface%par%alb_smax = alb_smax !0.79_dp !< max snow albedo
119 surface%par%alb_smin = alb_smin !0.6_dp !< min snow albedo
120 surface%par%albi = albi ! 0.41_dp !< albedo for ice
121 surface%par%albl = albl ! 0.07_dp !< albedo for land
122 surface%par%tmin = tmin ! -999_dp
123 surface%par%tmax = tmax ! 273.15_dp
124 surface%par%hcrit = hcrit !0.028_dp !< critical snow height for which grid cell is 50 % snow colvered [m]
125 surface%par%rcrit = rcrit !0.85_dp !< refreezing fraction is 50% [m]
126 surface%par%amp = Tamp !3.0_dp !< amplitude of diurnal cycle [K]
127 if (alb_scheme == 0) then
128 surface%par%alb_scheme="none"
129 else if (alb_scheme == 1) then
130 surface%par%alb_scheme = "slater"
131 else if (alb_scheme == 2) then
132 surface%par%alb_scheme = "denby"
133 else if (alb_scheme == 3) then
134 surface%par%alb_scheme = "isba"
135 else
136 print*, "ERROR: current albedo scheme is not available."
137 call exit(1)
138 end if
139 surface%par%tau_a = tau_a !0.008_dp
140 surface%par%tau_f = tau_f !0.24_dp
141 surface%par%w_crit = w_crit !15.0_dp ! default value
142 surface%par%mcrit = mcrit !6.0e-8_dp
143 surface%par%n_ksub = 3 ! sub ...
144 ! snow albedo of alex
145 surface%par%afac = afac
146 surface%par%tmid = tmid
147
148 ! initialize sub-daily time step tsticsub
149 surface%par%tsticsub = surface%par%tstic / dble(surface%par%n_ksub)
150
151 ! allocate necessary arrays for surface_physics module
152 call surface_alloc(surface%now,surface%par%nx)
153
154 ! initialise prognostic variables
155 if (debug) then
156 print*,"run_semic_transient: initialize variables."
157 end if
158 ! these values will be updated through "surface_energy_and_mass_balance" function.
159 surface%now%mask (:) = mask (:) ! 2.0_dp !loi_mask(:nx)
160 if (debug) then
161 print*,"run_semic_transient: initialize variables: mask"
162 end if
163 surface%now%hsnow (:) = hsnow (:) ! initial snow height...
164 surface%now%hice (:) = hice (:) ! initial ice height..
165 surface%now%tsurf (:) = tsurf_in (:) !< initial ice surface temperature
166 surface%now%alb (:) = albedo (:) !< initial albedo for energy balance.
167 surface%now%alb_snow(:) = albedo_snow(:) !< initial albedo for ISBA albedo method.
168 if (debug) then
169 print*,"run_semic_transient: initialize variables. DONE."
170 end if
171
172 if (debug) then
173 !print*, "====== global variable =========="
174 !print*, "nloop :", nloop
175 !print*, "nx :", surface%par%nx
176 !print*, "====== parameters ======"
177 !print*, "csh :", surface%par%csh
178 !print*, "clh :", surface%par%clh
179 !print*, "albeo scheme :", surface%par%alb_scheme
180 !print*, "albeo ice :", surface%par%albi
181 !print*, "tstic :", surface%par%tstic
182 !print*, "tsticsub :", surface%par%tsticsub
183 !print*, "n_ksub :", surface%par%n_ksub
184 print*, "====== inputs ========="
185 print*, "hsnow :", hsnow
186 print*, "====== state variables ======"
187 print*, "hsnow :", surface%now%hsnow
188 print*, "hice :", surface%now%hice
189 print*, "albeo :", surface%now%alb
190 print*, "albeo snow :", surface%now%alb_snow
191 print*, "mask :", surface%now%mask
192 print*, "tsurf :", surface%now%tsurf
193 print*, "sf :", sf_in
194 end if
195
196 ! define boundary conditions (not used, here!)
197 call surface_boundary_define(surface%bnd,surface%par%boundary)
198 !call print_boundary_opt(surface%bnd)
199
200 ! input with single value
201 do k =1,nloop
202 do i =1,ntime
203 if (debug) then
204 print*,"run_semic_transient: forcing data: ntime = ", i
205 end if
206 surface%now%sf = sf_in !(:,i)
207 surface%now%rf = rf_in !(:,i)
208 surface%now%sp = sp_in !(:,i)
209 surface%now%lwd = lwd_in !(:,i)
210 surface%now%swd = swd_in !(:,i)
211 surface%now%wind = wind_in!(:,i)
212 surface%now%rhoa = rhoa_in!(:,i)
213 surface%now%t2m = tt_in !(:,i)
214 surface%now%qq = qq_in !(:,i)
215 ! qmr_res is used to "energy_balance" in semic.
216 surface%now%qmr_res = qmr_in
217
218 ! calculate prognostic and diagnsotic variables
219 call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,i,year)
220
221 if (debug) then
222 print*,"done..."
223 end if
224 if (k == nloop) then
225 tsurf_out = surface%now%tsurf
226 ! melt - potential surface melt [m/s]
227 ! smb = SMB_ice + SMB_snow
228 ! smbi - SMB_ice (water equivalent m/sec)
229 ! smbs - SMB_snow (water equivalent m/sec)
230 smb_out =surface%now%smb ! smb = smb_snow + smb_ice
231 smbi_out =surface%now%smb_ice ! Csi (snow>ice) - melted_ice + refrezon_rain.
232 smbs_out =surface%now%smb_snow ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow
233 saccu_out =surface%now%acc ! acc = snowfall - sublimiation - refreezing
234 smelt_out =surface%now%melt ! potential surface melt = melt_ice + melt_snow
235 refr_out =surface%now%refr ! refreezing values. [m/sec]
236 alb_out =surface%now%alb
237 alb_snow_out =surface%now%alb_snow
238 hsnow_out =surface%now%hsnow
239 hice_out =surface%now%hice
240 qmr_out =surface%now%qmr_res
241 end if
242 end do
243 end do
244
245 ! de-allocate surface_physics arrays
246 call surface_dealloc(surface%now)
247
248end subroutine run_semic_transient ! }}}
Note: See TracBrowser for help on using the repository browser.