Skip to content
Snippets Groups Projects
Commit b42d39a0 authored by Victor Stepanenko's avatar Victor Stepanenko
Browse files

Dirichlet boundary condition for surface temperature implemented, not tested

parent 3a7df804
No related branches found
No related tags found
No related merge requests found
......@@ -67,6 +67,7 @@
real(kind=ireals), allocatable :: Short_net_rad(:)
real(kind=ireals), allocatable :: Sflux(:)
real(kind=ireals), allocatable :: cloud(:)
real(kind=ireals), allocatable :: SurfTemp(:)
real(kind=ireals), allocatable :: lat(:), lon(:)
......@@ -222,6 +223,7 @@
integer(kind=iintegers) :: N_SWdown
integer(kind=iintegers) :: N_LWdown
integer(kind=iintegers) :: N_NetRad
integer(kind=iintegers) :: N_SurfTemp
integer(kind=iintegers) :: N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud
integer(kind=iintegers) :: npoints, lakinterac, lakeform
integer(kind=iintegers) :: nstep_ncout
......@@ -278,7 +280,7 @@
& N_header_lines, N_coloumns, N_Year, N_Month, &
& N_Day, N_Hour, N_Precip, N_Uspeed, N_Vspeed, &
& N_Temp, N_Hum, N_Pres, N_SWdown, N_LWDown, N_NetRad, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, nstep_ncout, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, N_SurfTemp, nstep_ncout, &
& init_T, control_point, call_Flake, moving_average_window, mean_cycle_period, &
& nstep_out_Flake, dataname)
......@@ -316,6 +318,7 @@
allocate (Short_net_rad(1:npoints) )
allocate (Sflux(1:npoints) )
allocate (cloud(1:npoints) )
allocate (SurfTemp(1:npoints) )
allocate (lat(1:npoints))
allocate (lon(1:npoints))
......@@ -407,13 +410,15 @@
& kor, trib_inflow, effl_outflow, &
& widthchan, lengthchan, hbot, hbotchan, &
& Sals0, Salb0, fetch, phi, lam, us0, vs0, &
& Tm, alphax, alphay, a_veg, c_veg, h_veg, area_lake, cellipt, depth_area, &
& year0, month0, day0, npoints, lakinterac, lakeform, select_call, forc_format, form, &
& Tm, alphax, alphay, a_veg, c_veg, &
& h_veg, area_lake, cellipt, depth_area, &
& year0, month0, day0, npoints, lakinterac, &
& lakeform, select_call, forc_format, form, &
& spinup_times, spinup_period, cp_period, rad, &
& N_header_lines, N_coloumns, N_Year, N_Month, &
& N_Day, N_Hour, N_Precip, N_Uspeed, N_Vspeed, &
& N_Temp, N_Hum, N_Pres, N_SWdown, N_LWDown, N_NetRad, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, nstep_ncout, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, N_SurfTemp, nstep_ncout, &
& init_T, control_point, call_Flake, moving_average_window, mean_cycle_period, &
& nstep_out_Flake, dataname)
h10_2d(:,j) = h10
......@@ -595,9 +600,10 @@
& N_header_lines, N_coloumns, N_Year, N_Month, &
& N_Day, N_Hour, N_Precip, N_Uspeed, N_Vspeed, &
& N_Temp, N_Hum, N_Pres, N_SWdown, N_LWDown, N_NetRad, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, dataname, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, &
& N_cloud, N_SurfTemp, dataname, &
& ta, qa, pa, ua, va, atm_rad, Rad_bal, solar_rad, precip, &
& SensFlux,LatentFlux,Ustar, surfrad, cloud, &
& SensFlux, LatentFlux, Ustar, surfrad, cloud, SurfTemp, &
& outpath, spinup_done, npoints, dataread)
elseif (forc_format == 1) then ! Reading input file in netcdf format
......@@ -753,7 +759,7 @@
& alphay_2d(ix,iy), a_veg_2d(ix,iy), c_veg_2d(ix,iy), h_veg_2d(ix,iy), &
& area_lake_2d(ix,iy), cellipt_2d(ix,iy), depth_area_2d(1,1,ix,iy), &
& ts(ix),SensFlux(ix),LatentFlux(ix), exch_coef(ix), surfrad(ix), &
& cloud(ix), ftot(ix), h_2d(ix,iy), &
& cloud(ix), SurfTemp(ix), ftot(ix), h_2d(ix,iy), &
& ix,iy,1_iintegers,1_iintegers,nx,ny,nx,ny,ndatamax,year,month,day,hour, &
& init_T,flag_assim,flag_print, outpath, spinup_done, dataread, lakeform, &
& -999_iintegers,0_iintegers,(/0_iintegers,0_iintegers,0_iintegers/),.false., &
......@@ -1023,10 +1029,10 @@
& N_header_lines, N_coloumns, N_Year, N_Month, &
& N_Day, N_Hour, N_Precip, N_Uspeed, N_Vspeed, &
& N_Temp, N_Hum, N_Pres, N_SWdown, N_LWDown, N_NetRad, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, &
& N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud, N_SurfTemp, &
& dataname, &
& ta, qa, pa, ua, va, atm_rad, Rad_bal, solar_rad, precip, &
& SensFlux, LatentFlux, Ustar, surfrad, cloud, &
& SensFlux, LatentFlux, Ustar, surfrad, cloud, SurfTemp, &
& outpath,spinup_done,npoints,dataread)
use DRIVER_DATATYPES!, only : ireals, iintegers
......@@ -1074,6 +1080,7 @@
integer(kind=iintegers), intent(in) :: N_LWdown
integer(kind=iintegers), intent(in) :: N_NetRad
integer(kind=iintegers), intent(in) :: N_SensFlux, N_LatentFlux, N_Ustar, N_surfrad, N_cloud
integer(kind=iintegers), intent(in) :: N_SurfTemp
integer(kind=iintegers), intent(in) :: spinup_times
integer(kind=iintegers), intent(in) :: npoints
......@@ -1092,7 +1099,7 @@
real(kind=ireals), intent(out) :: solar_rad(1:npoints)
real(kind=ireals), intent(out) :: precip(1:npoints)
real(kind=ireals), intent(out) :: SensFlux(1:npoints), LatentFlux(1:npoints), &
& Ustar(1:npoints), surfrad(1:npoints), cloud(1:npoints)
& Ustar(1:npoints), surfrad(1:npoints), cloud(1:npoints), SurfTemp(1:npoints)
character(len=60), intent(out) :: outpath
......@@ -1115,6 +1122,7 @@
real(kind=ireals), allocatable :: Rad_balf(:)
real(kind=ireals), allocatable :: SensFluxf(:), LatentFluxf(:)
real(kind=ireals), allocatable :: Ustarf(:), surfradf(:), cloudf(:)
real(kind=ireals), allocatable :: SurfTempf(:)
real(kind=ireals), allocatable :: ta_old(:)
real(kind=ireals), allocatable :: qa_old(:)
......@@ -1127,6 +1135,7 @@
real(kind=ireals), allocatable :: precip_old(:)
real(kind=ireals), allocatable :: SensFlux_old(:), LatentFlux_old(:)
real(kind=ireals), allocatable :: Ustar_old(:), surfrad_old(:), cloud_old(:)
real(kind=ireals), allocatable :: SurfTemp_old(:)
real(kind=ireals) :: z0_wind
real(kind=ireals) :: sigma
......@@ -1188,6 +1197,7 @@
allocate(Ustarf(1:npoints))
allocate(surfradf(1:npoints))
allocate(cloudf(1:npoints))
allocate(SurfTempf(1:npoints))
allocate(ta_old(1:npoints))
allocate(qa_old(1:npoints))
......@@ -1203,6 +1213,7 @@
allocate(Ustar_old(1:npoints))
allocate(surfrad_old(1:npoints))
allocate(cloud_old(1:npoints))
allocate(SurfTemp_old(1:npoints))
if (form == 0) then
! This is a free form, adjusted in setup file
......@@ -1230,6 +1241,7 @@
Ustar_old (inpoints) = EXMISS(N_Ustar)
surfrad_old (inpoints) = EXMISS(N_surfrad)
cloud_old (inpoints) = EXMISS(N_cloud)
SurfTemp_old (inpoints) = EXMISS(N_SurfTemp)
!if (rad == 1) then
! atm_rad_old(inpoints) = xx
!elseif (rad == 2) then
......@@ -1255,6 +1267,7 @@
Ustarf(1:npoints) = Ustar_old(1:npoints)
surfradf(1:npoints) = surfrad_old(1:npoints)
cloudf(1:npoints) = cloud_old(1:npoints)
SurfTempf(1:npoints) = SurfTemp_old(1:npoints)
firstcallM1=.false.
......@@ -1290,6 +1303,7 @@
Ustar_old(:) = Ustarf(:)
surfrad_old(:) = surfradf(:)
cloud_old(:) = cloudf(:)
SurfTemp_old(:) = SurfTempf(:)
if (form == 0) then
do inpoints = 1, npoints
......@@ -1315,6 +1329,7 @@
Ustarf (inpoints) = EXMISS(N_Ustar)
surfradf (inpoints) = EXMISS(N_surfrad)
cloudf (inpoints) = EXMISS(N_cloud)
SurfTempf (inpoints) = EXMISS(N_SurfTemp)
!if (rad==1) then
! atm_radf(inpoints) = xx
!elseif (rad==2) then
......@@ -1359,6 +1374,8 @@
& (surfradf(inpoints) - surfrad_old(inpoints))/(hour_sec*interval)
cloud(inpoints) = cloud_old(inpoints) + (time - ti_old) * &
& (cloudf(inpoints) - cloud_old(inpoints))/(hour_sec*interval)
SurfTemp(inpoints) = SurfTemp_old(inpoints) + (time - ti_old) * &
& (SurfTempf(inpoints) - SurfTemp_old(inpoints))/(hour_sec*interval)
! The units of precipitation are m/sec
precip(inpoints) = precipf(inpoints)
......
......@@ -94,6 +94,7 @@
integer(kind=iintegers), private :: N_Ustar ; logical, private :: ok_N_Ustar = .false.
integer(kind=iintegers), private :: N_surfrad ; logical, private :: ok_N_surfrad = .false.
integer(kind=iintegers), private :: N_cloud ; logical, private :: ok_N_cloud = .false.
integer(kind=iintegers), private :: N_SurfTemp ; logical, private :: ok_N_SurfTemp = .false.
integer(kind=iintegers), private :: npoints ; logical, private :: ok_npoints = .false.
integer(kind=iintegers), private :: lakinterac ; logical, private :: ok_lakinterac = .false.
integer(kind=iintegers), private :: lakeform ; logical, private :: ok_lakeform = .false.
......@@ -129,7 +130,7 @@
& N_Day_outt, N_Hour_outt, N_Precip_outt, N_Uspeed_outt, N_Vspeed_outt, &
& N_Temp_outt, N_Hum_outt, N_Pres_outt, N_SWdown_outt, N_LWDown_outt, N_NetRad_outt, &
& N_SensFlux_outt, N_LatentFlux_outt, N_Ustar_outt, N_surfrad_outt, &
& N_cloud_outt,nstep_ncout_outt, &
& N_cloud_outt, N_SurfTemp_outt, nstep_ncout_outt, &
& init_T_outt, control_point_outt, call_Flake_outt, moving_average_window_outt, &
& mean_cycle_period_outt, nstep_out_Flake_outt, dataname_outt)
......@@ -211,6 +212,7 @@
integer(kind=iintegers), intent(out) :: N_Ustar_outt
integer(kind=iintegers), intent(out) :: N_surfrad_outt
integer(kind=iintegers), intent(out) :: N_cloud_outt
integer(kind=iintegers), intent(out) :: N_SurfTemp_outt
integer(kind=iintegers), intent(out) :: npoints_outt
integer(kind=iintegers), intent(out) :: lakinterac_outt
integer(kind=iintegers), intent(out) :: lakeform_outt
......@@ -398,6 +400,8 @@
& call WRITE_PAR_NOT_SET('N_surfrad')
if (.not.PAR_SET_INTEGER4(N_cloud, ok_N_cloud, N_cloud_outt)) &
& call WRITE_PAR_NOT_SET('N_cloud')
if (.not.PAR_SET_INTEGER4(N_SurfTemp, ok_N_SurfTemp, N_SurfTemp_outt)) &
& call WRITE_PAR_NOT_SET('N_SurfTemp')
if (.not.PAR_SET_INTEGER4(npoints, ok_npoints, npoints_outt)) &
& call WRITE_PAR_NOT_SET('npoints')
if (.not.PAR_SET_INTEGER4(lakinterac, ok_lakinterac, lakinterac_outt)) &
......@@ -714,6 +718,10 @@
N_cloud &
& = igetvarval(n1,n2,line,'N_cloud ')
ok_N_cloud = .true.
CASE ('N_SurfTemp')
N_SurfTemp &
& = igetvarval(n1,n2,line,'N_SurfTemp ')
ok_N_SurfTemp = .true.
CASE ('h10')
h10 = getvarval (n1,n2,line,'Init. depth of lake, m ')
ok_h10 = .true.
......
......@@ -146,7 +146,7 @@ SUBROUTINE LAKE &
& h10, l10, ls10, hs10, Ts0, Tb0, Tbb0, h_ML0, extwat, extice, &
& kor_in, trib_inflow, Sals0, Salb0, fetch, phi, lam, us0, vs0, &
& Tm, alphax, alphay, a_veg, c_veg, h_veg, area_lake, cellipt, depth_area, &
& tsw,hw1,xlew1,cdmw1,surfrad1,cloud1,ftot,h2_out, &
& tsw,hw1,xlew1,cdmw1,surfrad1,cloud1,SurfTemp1,ftot,h2_out, &
& ix,iy,nx0,ny0,nx,ny,nx_max,ny_max,ndatamax,year,month,day,hour,init_T,flag_assim,flag_print, &
& outpath1,spinup_done,dataread,lakeform,comm3d,rank_comm3d,coords,parallel,icp,step_final)
......@@ -304,6 +304,8 @@ real(kind=ireals), intent(in) :: precip1
real(kind=ireals), intent(in) :: Sflux01
!> Cloudiness, fraction
real(kind=ireals), intent(in) :: cloud1
!> Surface temperature, K
real(kind=ireals), intent(in) :: SurfTemp1
real(kind=ireals), intent(in) :: ch4_pres_atm, co2_pres_atm, o2_pres_atm
!> Time step for the LAKE model, s
......@@ -533,6 +535,7 @@ precip = precip1
Sflux0 = Sflux01*(roa0/row0)
outpath = outpath1
cloud = cloud1
SurfTemp = SurfTemp1
hw_input = hw1
xlew_input = xlew1
......@@ -929,7 +932,7 @@ if (multsoil) then
endif
! Calculation of the whole temperature profile
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
! Diagnostic calculations
do i = 1, M
......@@ -1313,7 +1316,7 @@ if (flag_snow == 1) then
endif
! Calculation of the whole temperature profile
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
!Assuming all dissolved species diffuse with the same molecular diffusivity as salinity
lamsal(:) = (lamw(:) - lamw0)/cw_m_row0 * alsal + lamsal0
......@@ -1580,7 +1583,7 @@ else
endif
! Calculation of the whole temperature profile
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
!Assuming all dissolved species diffuse with the same molecular diffusivity as salinity
lamsal(:) = (lamw(:) - lamw0)/cw_m_row0 * alsal + lamsal0
......@@ -1810,7 +1813,7 @@ if (flag_snow == 1) then
& gsp,soilflux(1,ix,iy),.false.,lsh)
! Temperature profile calculation
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
! Salinity profile
call S_DIFF(dt,ls,salice)
......@@ -1895,7 +1898,7 @@ else
& bathymwater,bathymice,bathymdice,bathymsoil, &
& gsp,soilflux(1,ix,iy),.false.,lsh)
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
! Salinity profile
call S_DIFF(dt,ls,salice)
......@@ -1951,7 +1954,7 @@ if4: IF (layer_case==4) THEN
call SOIL_COND_HEAT_COEF(nsoilcols)
endif
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
! Salinity profile
call S_DIFF(dt,ls,salice)
......@@ -1994,7 +1997,7 @@ if4: IF (layer_case==4) THEN
call SOIL_COND_HEAT_COEF(nsoilcols)
endif
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
! Salnity profile
call S_DIFF(dt,ls,salice)
......
......@@ -286,6 +286,7 @@ real(kind=ireals) :: velfrict, velfrict_prev, velfrict_bot
real(kind=ireals) :: cdm1
real(kind=ireals) :: turb_density_flux
real(kind=ireals) :: cloud
real(kind=ireals) :: SurfTemp
real(kind=ireals) :: discharge(1:2)
! Former common-block wind
......
......@@ -11,7 +11,7 @@ contains
!-----------------------------------------------------------------------------------
SUBROUTINE T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
& RadWater, RadIce, fetch, dt)
& RadWater, RadIce, SurfTemp, fetch, dt)
!T_SOLVER implements iterations to find water surface temperature
!at the next time step
......@@ -43,10 +43,12 @@ real(kind=ireals), intent(in) :: phi
!real(kind=ireals), intent(in) :: extwat, extice
real(kind=ireals), intent(in) :: fetch
real(kind=ireals), intent(in) :: dt
real(kind=ireals), intent(in) :: SurfTemp !> Surface temperature, Kelvin (prescribed, used as Dirichlet b.c.)
type(rad_type), intent(in) :: RadWater
type(rad_type), intent(in) :: RadIce
!Local variables
real(kind=ireals), allocatable :: Tsurf_prev(:,:)
integer(kind=iintegers), allocatable :: surftyp_prev(:,:)
......@@ -94,6 +96,15 @@ Tbc_if : if (cuette%par == 0) then
endif
call SURF_CHAR(surftyp(ix,iy),year,month,day,hour,phi,fetch)
Tbc_type : if (SurfTemp /= missing_value) then
!Using prescribed surface temperature as Dirichlet boundary condition
Tsurf3 = SurfTemp
call T_DIFF(1_iintegers,Tsurf3,dt,snow,ice,water,deepice,nsoilcols,.false.)
Bal3 = HEATBALANCE(Tsurf3,surftyp(ix,iy),RadWater,RadIce,fetch,dt)
else
if (nstep > 1.and.surftyp(ix,iy) == surftyp_prev(ix,iy)) then
Tsurf2 = Tsurf_prev(ix,iy) - 10.
!elseif (surftyp(ix,iy) == water_indic) then
......@@ -177,6 +188,7 @@ Tbc_if : if (cuette%par == 0) then
total_Bal = total_Bal + Bal3
endif Tbc_type
elseif (cuette%par == 1 .or. cuette%par == 11) then
! B.c.s for temperature for Cuette flow
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment