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