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