SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
anohm_module Module Reference

AnOHM: Analytical Objective Hysteresis Model. More...

Functions/Subroutines

subroutine anohm (tstep, dt_since_start, qn1, qn1_av_prev, dqndt_prev, qf, MetForcingData_grid, moist_surf, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, sfr, nsurf, EmissionsMethod, id, Gridiv, qn1_av_next, dqndt_next, a1, a2, a3, qs, deltaQi)
 High level wrapper for AnOHM calculation. More...
 
subroutine anohm_coef (sfc_typ, xid, xgrid, MetForcingData_grid, moist, EmissionsMethod, qf, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, xa1, xa2, xa3)
 High level wrapper for AnOHM coefficients calculation. More...
 
subroutine anohm_xts (sfc_typ, ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, xalb, xemis, xcp, xk, xch, xBo, tSd, xTHr, xTs)
 calculate the surface temperature related parameters (ATs, mTs, gamma) based on forcings and sfc. More...
 
subroutine anohm_coef_land_cal (ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, xalb, xemis, xcp, xk, xch, xBo, xa1, xa2, xa3, ATs, mTs, gamma)
 
subroutine anohm_coef_water_cal (ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, xalb, xemis, xcp, xk, xch, xBo, xeta, xmu, xa1, xa2, xa3, ATs, mTs, gamma)
 a wrapper for retrieving AnOHM coefficients of water body More...
 
subroutine anohm_fc (xid, MetForcingData_grid, EmissionsMethod, qf, ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
 
subroutine anohm_fcload (xid, MetForcingData_grid, EmissionsMethod, qf, Sd, Ta, RH, pres, WS, WF, AH, tHr)
 load forcing series for AnOHM_FcCal More...
 
subroutine anohm_fccal (Sd, Ta, WS, WF, AH, tHr, ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
 calculate the key parameters of a sinusoidal curve for AnOHM forcings i.e., a, b, c in a*Sin(Pi/12*t+b)+c More...
 
subroutine anohm_shapefit (tHr, obs, amp, mean, tpeak)
 calculate the key parameters of a sinusoidal curve for AnOHM forcings i.e., a, b, c in a*Sin(Pi/12*t+b)+c, where t is in hour More...
 
subroutine fsin (m, n, x, xdat, ydat, fvec, iflag)
 sinusoidal function f(t) for fitting: f(t) = mean+amp*Sin(Pi/12(t-delta)) x = (/mean,amp,delta/) contains the fitting parameters More...
 
subroutine anohm_bo_cal (sfc_typ, Sd, Ta, RH, pres, tHr, ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, xalb, xemis, xcp, xk, xch, xSM, tSd, xBo)
 estimate daytime Bowen ratio for calculation of AnOHM coefficients More...
 
subroutine fcnbo (n, x, fvec, iflag, m, prms)
 this fucntion will construct an equaiton for Bo calculation More...
 
real(kind(1d0)) function esat_fn (Ta)
 calculate saturation vapor pressure (es) at air temperature (Ta) (MRR, 1987) More...
 
real(kind(1d0)) function qsat_fn (Ta, pres)
 calculate saturation specific humidity (qsat) at air temperature (Ta) and atmospheric pressure (pres) (MRR, 1987) More...
 
real(kind(1d0)) function qa_fn (Ta, RH, pres)
 convert relative humidity (RH) to specific humidity (qa) at air temperature (Ta) and atmospheric pressure (pres) More...
 

Detailed Description

AnOHM: Analytical Objective Hysteresis Model.

Author
Ting Sun, ting..nosp@m.sun@.nosp@m.readi.nosp@m.ng.a.nosp@m.c.uk

calculate heat storage. model details refer to https://doi.org/10.5194/gmd-2016-300

Function/Subroutine Documentation

◆ anohm()

subroutine anohm_module::anohm ( integer, intent(in)  tstep,
integer, intent(in)  dt_since_start,
real(kind(1d0)), intent(in)  qn1,
real(kind(1d0)), intent(in)  qn1_av_prev,
real(kind(1d0)), intent(in)  dqndt_prev,
real(kind(1d0)), intent(in)  qf,
real(kind(1d0)), dimension(:, :), intent(in)  MetForcingData_grid,
real(kind(1d0)), dimension(nsurf), intent(in)  moist_surf,
real(kind(1d0)), dimension(:), intent(in)  alb,
real(kind(1d0)), dimension(:), intent(in)  emis,
real(kind(1d0)), dimension(:), intent(in)  cpAnOHM,
real(kind(1d0)), dimension(:), intent(in)  kkAnOHM,
real(kind(1d0)), dimension(:), intent(in)  chAnOHM,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
integer, intent(in)  nsurf,
integer, intent(in)  EmissionsMethod,
integer, intent(in)  id,
integer, intent(in)  Gridiv,
real(kind(1d0)), intent(out)  qn1_av_next,
real(kind(1d0)), intent(out)  dqndt_next,
real(kind(1d0)), intent(out)  a1,
real(kind(1d0)), intent(out)  a2,
real(kind(1d0)), intent(out)  a3,
real(kind(1d0)), intent(out)  qs,
real(kind(1d0)), dimension(nsurf), intent(out)  deltaQi 
)

High level wrapper for AnOHM calculation.

calculate heat storage based within AnOHM framework.

Returns
  1. grid ensemble heat storage: QS = a1*(Q*)+a2*(dQ*/dt)+a3
  2. grid ensemble OHM coefficients: a1, a2 and a3
Parameters
[in]metforcingdata_gridmet forcing array of grid
[in]qn1net all-wave radiation [W m-2]
[in]qfanthropogenic heat flux [W m-2]
[in]sfrsurface fraction (0-1) [-]
[in]moist_surfnon-dimensional surface wetness status (0-1) [-]
[in]albalbedo [-]
[in]emisemissivity [-]
[in]cpanohmheat capacity [J m-3 K-1]
[in]kkanohmthermal conductivity [W m-1 K-1]
[in]chanohmbulk transfer coef [J m-3 K-1]
[in]idday of year [-]
[in]gridivgrid id [-]
[in]emissionsmethodAnthropHeat option [-]
[in]nsurfnumber of surfaces [-]
[out]a1AnOHM coefficients of grid [-]
[out]a2AnOHM coefficients of grid [h]
[out]a3AnOHM coefficients of grid [W m-2]
[out]qsstorage heat flux [W m-2]
[out]deltaqistorage heat flux of snow surfaces

Definition at line 36 of file suews_phys_anohm.f95.

References anohm_coef(), errorhint(), ohm_dqndt_cal_x(), and ohm_qs_cal().

Referenced by suews_driver::suews_cal_qs().

36 
37  IMPLICIT NONE
38  INTEGER, INTENT(in) :: tstep ! time step [s]
39  INTEGER, INTENT(in) :: dt_since_start ! time since simulation starts [s]
40 
41  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :)::metforcingdata_grid
42 
43  REAL(KIND(1d0)), INTENT(in):: qn1
44  REAL(KIND(1d0)), INTENT(in):: qf
45  REAL(KIND(1d0)), INTENT(in):: sfr(nsurf)
46  REAL(KIND(1d0)), INTENT(in):: moist_surf(nsurf)
47 
48  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::alb
49  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::emis
50  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::cpanohm
51  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::kkanohm
52  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::chanohm
53 
54  INTEGER, INTENT(in):: id
55  INTEGER, INTENT(in):: gridiv
56  INTEGER, INTENT(in):: emissionsmethod
57  INTEGER, INTENT(in):: nsurf
58  ! INTEGER,INTENT(in):: nsh !< number of timesteps in one hour [-]
59 
60  REAL(KIND(1d0)), INTENT(in)::qn1_av_prev
61  REAL(KIND(1d0)), INTENT(out)::qn1_av_next
62  REAL(KIND(1d0)), INTENT(in)::dqndt_prev !Rate of change of net radiation [W m-2 h-1] at t-1
63  REAL(KIND(1d0)), INTENT(out)::dqndt_next !Rate of change of net radiation [W m-2 h-1] at t-1
64 
65  ! REAL(KIND(1d0)),INTENT(inout)::qn1_store(nsh) !< stored qn1 [W m-2]
66  ! REAL(KIND(1d0)),INTENT(inout)::qn1_av_store(2*nsh+1) !< average net radiation over previous hour [W m-2]
67 
68  REAL(KIND(1d0)), INTENT(out):: a1
69  REAL(KIND(1d0)), INTENT(out):: a2
70  REAL(KIND(1d0)), INTENT(out):: a3
71  REAL(KIND(1d0)), INTENT(out):: qs
72  REAL(KIND(1d0)), INTENT(out):: deltaqi(nsurf)
73 
74  INTEGER :: is, xid
75  INTEGER, SAVE :: id_save ! store index of the valid day with enough data
76  REAL(KIND(1d0)), PARAMETER::notused = -55.5
77  INTEGER, PARAMETER::notusedi = -55
78  LOGICAL :: idq ! whether id contains enough data
79 
80  ! REAL(KIND(1d0)) :: dqndt !< rate of change of net radiation [W m-2 h-1] at t-2
81  ! REAL(KIND(1d0)) :: surfrac !< surface fraction accounting for SnowFrac if appropriate
82  REAL(KIND(1d0)), DIMENSION(nsurf) :: xa1, xa2, xa3
83  ! REAL(KIND(1d0)) :: qn1_av ! average net radiation over previous hour [W m-2]
84  ! REAL(KIND(1d0)) :: nsh_nna ! number of timesteps per hour with non -999 values (used for spinup)
85 
86  ! initialize output variables
87  xa1 = 0.1
88  xa2 = 0.2
89  xa3 = 10
90  qs = -999
91  deltaqi = 0 ! NB: as snow part is not implemented within AnOHM yet, this is just a placeholder
92 
93  ! to test if the current met block contains enough data for AnOHM
94  ! TODO: more robust selection should be implemented
95  ! daylight hours >= 6
96  idq = count(metforcingdata_grid(:, 2) == id .AND. & ! day of year
97  metforcingdata_grid(:, 4) == 0 .AND. & ! minutes
98  metforcingdata_grid(:, 15) > 0) & ! Sd
99  >= 6
100 
101  ! PRINT*, idQ
102  IF (idq) THEN
103  ! given enough data, calculate coefficients of day `id`
104  xid = id
105  id_save = id ! store index of the valid day with enough data
106  ELSE
107  ! otherwise calculate coefficients of yesterday: `id-1`
108  xid = id_save
109  END IF
110 
111  DO is = 1, nsurf
112  ! call AnOHM to calculate the coefs.
113  CALL anohm_coef(is, xid, gridiv, metforcingdata_grid, moist_surf, emissionsmethod, qf, & !input
114  alb, emis, cpanohm, kkanohm, chanohm, &! input
115  xa1(is), xa2(is), xa3(is)) ! output
116  ! print*, 'AnOHM_coef are: ',xa1,xa2,xa3
117  ENDDO
118 
119  ! calculate the areally-weighted OHM coefficients
120  a1 = dot_product(xa1, sfr)
121  a2 = dot_product(xa2, sfr)
122  a3 = dot_product(xa3, sfr)
123 
124  ! Calculate radiation part ------------------------------------------------------------
125  qs = -999 !qs = Net storage heat flux [W m-2]
126  IF (qn1 > -999) THEN !qn1 = Net all-wave radiation [W m-2]
127 
128  ! Store instantaneous qn1 values for previous hour (qn1_store) and average (qn1_av)
129  ! CALL OHM_dqndt_cal(nsh,qn1,qn1_store,qn1_av_store,dqndt)
130  CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn1_av_prev, qn1, dqndt_prev, &
131  qn1_av_next, dqndt_next)
132 
133  ! Calculate net storage heat flux
134  CALL ohm_qs_cal(qn1, dqndt_prev, a1, a2, a3, qs)
135 
136  ELSE
137  CALL errorhint(21, 'SUEWS_AnOHM.f95: bad value for qn found during qs calculation.', qn1, notused, notusedi)
138  ENDIF
139 
allocatearray::metforcingdata_grid
real(kind(1d0)), dimension(:, :), allocatable metforcingdata_grid
Definition: suews_ctrl_const.f95:116
allocatearray::a2
real(kind(1d0)) a2
Definition: suews_ctrl_const.f95:394
defaultnotused::notused
real(kind(1d0)) notused
Definition: suews_ctrl_const.f95:1244
allocatearray::a1
real(kind(1d0)) a1
Definition: suews_ctrl_const.f95:394
allocatearray::kkanohm
real(kind(1d0)), dimension(nsurf) kkanohm
Definition: suews_ctrl_const.f95:460
allocatearray::nsurf
integer, parameter nsurf
Definition: suews_ctrl_const.f95:183
allocatearray::a3
real(kind(1d0)) a3
Definition: suews_ctrl_const.f95:394
data_in::emissionsmethod
integer emissionsmethod
Definition: suews_ctrl_const.f95:977
data_in::qn1
real(kind(1d0)) qn1
Definition: suews_ctrl_const.f95:1042
time::id
integer id
Definition: suews_ctrl_const.f95:1252
allocatearray::chanohm
real(kind(1d0)), dimension(nsurf) chanohm
Definition: suews_ctrl_const.f95:461
allocatearray::emis
real(kind(1d0)), dimension(nsurf) emis
Definition: suews_ctrl_const.f95:361
sues_data::tstep
integer tstep
Definition: suews_ctrl_const.f95:1384
time::dt_since_start
integer dt_since_start
Definition: suews_ctrl_const.f95:1262
allocatearray::deltaqi
real(kind(1d0)), dimension(nsurf) deltaqi
Definition: suews_ctrl_const.f95:410
sues_data::is
integer is
Definition: suews_ctrl_const.f95:1403
allocatearray::cpanohm
real(kind(1d0)), dimension(nsurf) cpanohm
Definition: suews_ctrl_const.f95:459
allocatearray::sfr
real(kind(1d0)), dimension(nsurf) sfr
Definition: suews_ctrl_const.f95:200
ohm_dqndt_cal_x
subroutine ohm_dqndt_cal_x(dt, dt_since_start, qn1_av_prev, qn1, dqndt_prev, qn1_av_next, dqndt_next)
Definition: suews_phys_ohm.f95:272
allocatearray::alb
real(kind(1d0)), dimension(nsurf) alb
Definition: suews_ctrl_const.f95:360
errorhint
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
Definition: suews_ctrl_error.f95:2
ohm_qs_cal
subroutine ohm_qs_cal(qn1, dqndt, a1, a2, a3, qs)
Definition: suews_phys_ohm.f95:348
Here is the call graph for this function:
Here is the caller graph for this function:

◆ anohm_bo_cal()

subroutine anohm_module::anohm_bo_cal ( integer, intent(in)  sfc_typ,
real(kind=8), dimension(:), intent(in)  Sd,
real(kind=8), dimension(:), intent(in)  Ta,
real(kind=8), dimension(:), intent(in)  RH,
real(kind=8), dimension(:), intent(in)  pres,
real(kind=8), dimension(:), intent(in)  tHr,
real(kind(1d0)), intent(in)  ASd,
real(kind(1d0)), intent(in)  mSd,
real(kind(1d0)), intent(in)  ATa,
real(kind(1d0)), intent(in)  mTa,
real(kind(1d0)), intent(in)  tau,
real(kind(1d0)), intent(in)  mWS,
real(kind(1d0)), intent(in)  mWF,
real(kind(1d0)), intent(in)  mAH,
real(kind(1d0)), intent(in)  xalb,
real(kind(1d0)), intent(in)  xemis,
real(kind(1d0)), intent(in)  xcp,
real(kind(1d0)), intent(in)  xk,
real(kind(1d0)), intent(in)  xch,
real(kind(1d0)), intent(in)  xSM,
real(kind(1d0)), intent(in)  tSd,
real(kind=8), intent(out)  xBo 
)

estimate daytime Bowen ratio for calculation of AnOHM coefficients

Parameters
[in]sdincoming solar radiation [W m-2]
[in]taair temperature [degC]
[in]rhrelative humidity [%]
[in]presAtmospheric pressure [mbar]
[in]thrlocal time [hr]
[in]asddaily amplitude of solar radiation [W m-2]
[in]msddaily mean solar radiation [W m-2]
[in]tsdlocal peaking time of solar radiation [hr]
[in]atadaily amplitude of air temperature [degC]
[in]mtadaily mean air temperature [degC]
[in]tauphase lag between Sd and Ta (Ta-Sd) [rad]
[in]mwsdaily mean wind speed [m s-1]
[in]mwfdaily mean underground moisture flux [m3 s-1 m-2]
[in]mahdaily mean anthropogenic heat flux [W m-2]
[in]xalbalbedo [-]
[in]xemisemissivity [-]
[in]xcpheat capacity [J m-3 K-1]
[in]xkthermal conductivity [W m-1 K-1]
[in]xchbulk transfer coef [J m-3 K-1]
[in]xsmsurface moisture status [-]

Definition at line 1044 of file suews_phys_anohm.f95.

References fcnbo(), and hybrd1().

Referenced by anohm_coef().

1044 
1045  IMPLICIT NONE
1046 
1047  ! input:
1048  INTEGER, INTENT(in) :: sfc_typ ! unknown Bowen ratio
1049 
1050  ! input: daytime series
1051  REAL(kind=8), INTENT(in), DIMENSION(:) ::sd
1052  REAL(kind=8), INTENT(in), DIMENSION(:) ::ta
1053  REAL(kind=8), INTENT(in), DIMENSION(:) ::rh
1054  REAL(kind=8), INTENT(in), DIMENSION(:) ::pres
1055  REAL(kind=8), INTENT(in), DIMENSION(:) ::thr
1056 
1057  ! input: forcing scales
1058  REAL(KIND(1d0)), INTENT(in):: asd
1059  REAL(KIND(1d0)), INTENT(in):: msd
1060  REAL(KIND(1d0)), INTENT(in):: tsd
1061  REAL(KIND(1d0)), INTENT(in):: ata
1062  REAL(KIND(1d0)), INTENT(in):: mta
1063  REAL(KIND(1d0)), INTENT(in):: tau
1064  REAL(KIND(1d0)), INTENT(in):: mws
1065  REAL(KIND(1d0)), INTENT(in):: mwf
1066  REAL(KIND(1d0)), INTENT(in):: mah
1067 
1068  ! input: surface properties
1069  REAL(KIND(1d0)), INTENT(in) :: xalb
1070  REAL(KIND(1d0)), INTENT(in) :: xemis
1071  REAL(KIND(1d0)), INTENT(in) :: xcp
1072  REAL(KIND(1d0)), INTENT(in) :: xk
1073  REAL(KIND(1d0)), INTENT(in) :: xch
1074  REAL(KIND(1d0)), INTENT(in) :: xsm
1075 
1076  ! output:
1077  REAL(kind=8), INTENT(out) :: xbo ! unknown Bowen ratio [-]
1078 
1079  ! EXTERNAL fcnBo
1080 
1081  REAL(kind=8), ALLOCATABLE :: x(:), fvec(:), prms(:)
1082 
1083  INTEGER :: lenday, n, m, info, err, nvar, nprm
1084  LOGICAL, DIMENSION(:), ALLOCATABLE :: maskday
1085  REAL(kind=8) :: tol = 1e-20
1086 
1087  ! LOGICAL, ALLOCATABLE,DIMENSION(:) :: metMask
1088 
1089  ! daytime mask
1090  ALLOCATE (maskday(SIZE(sd)), stat=err)
1091  maskday = sd > 5
1092  ! length of daytime series
1093  lenday = SIZE(pack(sd, mask=maskday), dim=1)
1094  ! ! daytime mask
1095  ! IF (ALLOCATED(metMask)) DEALLOCATE(metMask, stat=err)
1096  ! ALLOCATE(metMask(lenDay))
1097  ! metMask=(Sd>0)
1098 
1099  ! assign x vector:
1100  n = 1
1101  ALLOCATE (x(n), stat=err)
1102  IF (err /= 0) print *, "x: Allocation request denied"
1103  ALLOCATE (fvec(n), stat=err)
1104  IF (err /= 0) print *, "fvec: Allocation request denied"
1105 
1106  ! pass initial Bowen ratio:
1107  xbo = 10
1108  x(1) = xbo
1109 
1110  ! NB: set these numbers properly if any changes made to this subroutine:
1111  ! number of parameters to pass:
1112  nprm = 16
1113  ! number of met variables to pass:
1114  nvar = 5
1115  ! length of x vector that holds unknown and parameters
1116  m = nprm + nvar*lenday
1117 
1118  ALLOCATE (prms(m), stat=err)
1119  IF (err /= 0) print *, "prms: Allocation request denied"
1120 
1121  ! pass forcing scales:
1122  prms(1) = asd
1123  prms(2) = msd
1124  prms(3) = ata
1125  prms(4) = mta
1126  prms(5) = tau
1127  prms(6) = mws
1128  prms(7) = mwf
1129  prms(8) = mah
1130 
1131  ! pass sfc. property scales:
1132  prms(9) = xalb
1133  prms(10) = xemis
1134  prms(11) = xcp
1135  prms(12) = xk
1136  prms(13) = xch
1137  prms(14) = xsm
1138 
1139  ! pass tSd:
1140  prms(15) = tsd
1141 
1142  ! pass tSd:
1143  prms(16) = sfc_typ*1.0
1144 
1145  ! extract daytime series
1146  prms(nprm + 1:m) = pack((/sd, ta, rh, pres, thr/), &
1147  mask=pack(spread(maskday, dim=2, ncopies=nvar), .true.))
1148 
1149  ! PRINT*, 'xBo before solve:',x(1)
1150  ! PRINT*, 'fvec before solve:',fvec(1)
1151  ! PRINT*, 'xSM:',xSM
1152  ! solve nonlinear equation fcnBo(x)=0
1153  CALL hybrd1(fcnbo, n, x, fvec, tol, info, m, prms)
1154  xbo = x(1)
1155  ! PRINT*, 'xBo after solve: ',x(1)
1156  ! PRINT*, 'fvec after solve:',fvec(1)
1157 
1158  IF (ALLOCATED(x)) DEALLOCATE (x, stat=err)
1159  IF (err /= 0) print *, "x: Deallocation request denied"
1160  IF (ALLOCATED(fvec)) DEALLOCATE (fvec, stat=err)
1161  IF (err /= 0) print *, "fvec: Deallocation request denied"
1162  IF (ALLOCATED(prms)) DEALLOCATE (prms, stat=err)
1163  IF (err /= 0) print *, "prms: Deallocation request denied"
1164 
hybrd1
subroutine hybrd1(fcn, n, x, fvec, tol, info, m, prms)
Definition: suews_util_minpack.f95:1414
allocatearray::xbo
real(kind(1d0)) xbo
Definition: suews_ctrl_const.f95:320
Here is the call graph for this function:
Here is the caller graph for this function:

◆ anohm_coef()

subroutine anohm_module::anohm_coef ( integer, intent(in)  sfc_typ,
integer, intent(in)  xid,
integer, intent(in)  xgrid,
real(kind(1d0)), dimension(:, :), intent(in)  MetForcingData_grid,
real(kind(1d0)), dimension(:), intent(in)  moist,
integer, intent(in)  EmissionsMethod,
real(kind(1d0)), intent(in)  qf,
real(kind(1d0)), dimension(:), intent(in)  alb,
real(kind(1d0)), dimension(:), intent(in)  emis,
real(kind(1d0)), dimension(:), intent(in)  cpAnOHM,
real(kind(1d0)), dimension(:), intent(in)  kkAnOHM,
real(kind(1d0)), dimension(:), intent(in)  chAnOHM,
real(kind(1d0)), intent(out)  xa1,
real(kind(1d0)), intent(out)  xa2,
real(kind(1d0)), intent(out)  xa3 
)

High level wrapper for AnOHM coefficients calculation.

calculate OHM coefficients within AnOHM framework.

Returns
  1. OHM coefficients of a given surface type: a1, a2 and a3
Parameters
[in]sfc_typsurface type [-]
[in]xidday of year [-]
[in]xgridgrid id [-]
[in]emissionsmethodAnthropHeat option [-]
[in]qfanthropogenic heat flux [W m-2]
[in]albalbedo [-]
[in]emisemissivity [-]
[in]cpanohmheat capacity [J m-3 K-1]
[in]kkanohmthermal conductivity [W m-1 K-1]
[in]chanohmbulk transfer coef [J m-3 K-1]
[in]moistsurface wetness status [-]
[in]metforcingdata_gridmet forcing array of grid
[out]xa1AnOHM coefficients of grid [-]
[out]xa2AnOHM coefficients of grid [h]
[out]xa3AnOHM coefficients of grid [W m-2]

Definition at line 154 of file suews_phys_anohm.f95.

References anohm_bo_cal(), anohm_coef_land_cal(), anohm_coef_water_cal(), anohm_fc(), and anohm_fcload().

Referenced by anohm().

154 
155  IMPLICIT NONE
156 
157  ! input
158  INTEGER, INTENT(in):: sfc_typ
159  INTEGER, INTENT(in):: xid
160  INTEGER, INTENT(in):: xgrid
161  INTEGER, INTENT(in):: emissionsmethod
162 
163  REAL(KIND(1d0)), INTENT(in):: qf
164  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: alb
165  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: emis
166  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: cpanohm
167  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: kkanohm
168  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: chanohm
169  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: moist
170  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :) :: metforcingdata_grid
171 
172  ! output
173  REAL(KIND(1d0)), INTENT(out) :: xa1
174  REAL(KIND(1d0)), INTENT(out) :: xa2
175  REAL(KIND(1d0)), INTENT(out) :: xa3
176 
177  ! surface temperature related scales:
178  REAL(KIND(1d0)):: ats
179  REAL(KIND(1d0)):: mts
180  REAL(KIND(1d0)):: gamma
181 
182  ! forcing scales
183  REAL(KIND(1d0))::asd, msd, tsd
184  REAL(KIND(1d0))::ata, mta, tta
185  REAL(KIND(1d0))::tau
186  REAL(KIND(1d0))::mws, mwf, mah
187 
188  ! forcings:
189  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::sd ! incoming solar radiation [W m-2]
190  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::ta ! air temperature [degC]
191  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::rh ! relative humidity [%]
192  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::pres ! air pressure [hPa]
193  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::ws ! wind speed [m s-1]
194  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::wf ! water flux density [m3 s-1 m-2]
195  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::ah ! anthropogenic heat [W m-2]
196  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::thr ! time [hr]
197 
198  ! sfc. properties:
199  REAL(KIND(1d0)) ::xalb ! albedo,
200  REAL(KIND(1d0)) ::xemis ! emissivity,
201  REAL(KIND(1d0)) ::xcp ! heat capacity,
202  REAL(KIND(1d0)) ::xk ! thermal conductivity,
203  REAL(KIND(1d0)) ::xch ! bulk transfer coef.
204  REAL(KIND(1d0)) ::xbo ! Bowen ratio
205  REAL(KIND(1d0)) ::xeta ! effective absorption coefficient
206  REAL(KIND(1d0)) ::xmu ! effective absorption fraction
207  REAL(KIND(1d0)) ::xmoist ! surface wetness
208 
209  ! locally saved variables:
210  ! if coefficients have been calculated, just reload them
211  ! otherwise, do the calculation
212  INTEGER, SAVE :: id_save, grid_save
213  REAL(KIND(1d0)), SAVE:: coeff_grid_day(7, 3) = -999.
214 
215  ! PRINT*, 'xid,id_save',xid,id_save
216  ! PRINT*, 'xgrid,grid_save',xgrid,grid_save
217  ! PRINT*, 'sfc_typ',sfc_typ
218  ! PRINT*, 'coeff_grid_day',coeff_grid_day(sfc_typ,:)
219  IF (xid == id_save .AND. xgrid == grid_save) THEN
220  ! if coefficients have been calculated, just reload them
221  ! print*, 'here no repetition'
222  xa1 = coeff_grid_day(sfc_typ, 1)
223  xa2 = coeff_grid_day(sfc_typ, 2)
224  xa3 = coeff_grid_day(sfc_typ, 3)
225  ELSE
226 
227  ! load forcing characteristics:
228  CALL anohm_fc( &
229  xid, metforcingdata_grid, emissionsmethod, qf, & ! input
230  asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah) ! output
231 
232  ! load forcing variables:
233  CALL anohm_fcload( &
234  xid, metforcingdata_grid, emissionsmethod, qf, & ! input
235  sd, ta, rh, pres, ws, wf, ah, thr) ! output
236 
237  ! load sfc. properties:
238  xalb = alb(sfc_typ)
239  xemis = emis(sfc_typ)
240  xcp = cpanohm(sfc_typ)
241  xk = kkanohm(sfc_typ)
242  xch = chanohm(sfc_typ)
243  xmoist = moist(sfc_typ)
244 
245  ! PRINT*, 'xBo before:',xBo
246  ! calculate Bowen ratio:
247  CALL anohm_bo_cal( &
248  sfc_typ, &
249  sd, ta, rh, pres, thr, & ! input: forcing
250  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
251  xalb, xemis, xcp, xk, xch, xmoist, & ! input: sfc properties
252  tsd, & ! input: peaking time of Sd in hour
253  xbo) ! output: Bowen ratio
254  ! PRINT*, 'xBo after:',xBo
255 
256  ! calculate AnOHM coefficients
257  SELECT CASE (sfc_typ)
258  CASE (1:6) ! land surfaces
259  CALL anohm_coef_land_cal( &
260  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
261  xalb, xemis, xcp, xk, xch, xbo, & ! input: sfc properties
262  xa1, xa2, xa3, ats, mts, gamma) ! output: surface temperature related scales by AnOHM
263 
264  CASE (7) ! water surface
265  ! NB:give fixed values for the moment
266  xeta = 0.3
267  xmu = 0.2
268  CALL anohm_coef_water_cal( &
269  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
270  xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, & ! input: sfc properties
271  xa1, xa2, xa3, ats, mts, gamma) ! output
272 
273  ! save variables for marking status as done
274  id_save = xid
275  grid_save = xgrid
276 
277  END SELECT
278  ! save variables for reusing values of the same day
279  coeff_grid_day(sfc_typ, :) = (/xa1, xa2, xa3/)
280  END IF
281 
allocatearray::metforcingdata_grid
real(kind(1d0)), dimension(:, :), allocatable metforcingdata_grid
Definition: suews_ctrl_const.f95:116
estm_data::ws
real(kind(1d0)) ws
Definition: suews_phys_estm.f95:85
allocatearray::kkanohm
real(kind(1d0)), dimension(nsurf) kkanohm
Definition: suews_ctrl_const.f95:460
data_in::emissionsmethod
integer emissionsmethod
Definition: suews_ctrl_const.f95:977
allocatearray::chanohm
real(kind(1d0)), dimension(nsurf) chanohm
Definition: suews_ctrl_const.f95:461
allocatearray::emis
real(kind(1d0)), dimension(nsurf) emis
Definition: suews_ctrl_const.f95:361
allocatearray::cpanohm
real(kind(1d0)), dimension(nsurf) cpanohm
Definition: suews_ctrl_const.f95:459
allocatearray::alb
real(kind(1d0)), dimension(nsurf) alb
Definition: suews_ctrl_const.f95:360
allocatearray::xbo
real(kind(1d0)) xbo
Definition: suews_ctrl_const.f95:320
Here is the call graph for this function:
Here is the caller graph for this function:

◆ anohm_coef_land_cal()

subroutine anohm_module::anohm_coef_land_cal ( real(kind(1d0)), intent(in)  ASd,
real(kind(1d0)), intent(in)  mSd,
real(kind(1d0)), intent(in)  ATa,
real(kind(1d0)), intent(in)  mTa,
real(kind(1d0)), intent(in)  tau,
real(kind(1d0)), intent(in)  mWS,
real(kind(1d0)), intent(in)  mWF,
real(kind(1d0)), intent(in)  mAH,
real(kind(1d0)), intent(in)  xalb,
real(kind(1d0)), intent(in)  xemis,
real(kind(1d0)), intent(in)  xcp,
real(kind(1d0)), intent(in)  xk,
real(kind(1d0)), intent(in)  xch,
real(kind(1d0)), intent(in)  xBo,
real(kind(1d0)), intent(out)  xa1,
real(kind(1d0)), intent(out)  xa2,
real(kind(1d0)), intent(out)  xa3,
real(kind(1d0)), intent(out)  ATs,
real(kind(1d0)), intent(out)  mTs,
real(kind(1d0)), intent(out)  gamma 
)
Parameters
[in]asddaily amplitude of solar radiation [W m-2]
[in]msddaily mean solar radiation [W m-2]
[in]atadaily amplitude of air temperature [K]
[in]mtadaily mean air temperature [K]
[in]tauphase lag between Sd and Ta (Ta-Sd) [rad]
[in]mwsdaily mean wind speed [m s-1]
[in]mwfdaily mean underground moisture flux [m3 s-1 m-2]
[in]mahdaily mean anthropogenic heat flux [W m-2]
[in]xalbalbedo [-]
[in]xemisemissivity [-]
[in]xcpheat capacity [J m-3 K-1]
[in]xkthermal conductivity [W m-1 K-1]
[in]xchbulk transfer coef [J m-3 K-1]
[in]xboBowen ratio [-]
[out]xa1AnOHM coefficients of grid [-]
[out]xa2AnOHM coefficients of grid [h]
[out]xa3AnOHM coefficients of grid [W m-2]
[out]atsdaily amplitude of surface temperature [K]
[out]mtsdaily mean of surface temperature [K]
[out]gammaphase difference between Ts and Sd [K]

Definition at line 367 of file suews_phys_anohm.f95.

Referenced by anohm_coef(), and anohm_xts().

367 
368  IMPLICIT NONE
369 
370  ! input: forcing scales
371  REAL(KIND(1d0)), INTENT(in):: asd
372  REAL(KIND(1d0)), INTENT(in):: msd
373  REAL(KIND(1d0)), INTENT(in):: ata
374  REAL(KIND(1d0)), INTENT(in):: mta
375  REAL(KIND(1d0)), INTENT(in):: tau
376  REAL(KIND(1d0)), INTENT(in):: mws
377  REAL(KIND(1d0)), INTENT(in):: mwf
378  REAL(KIND(1d0)), INTENT(in):: mah
379  ! input: sfc properties
380  REAL(KIND(1d0)), INTENT(in):: xalb
381  REAL(KIND(1d0)), INTENT(in):: xemis
382  REAL(KIND(1d0)), INTENT(in):: xcp
383  REAL(KIND(1d0)), INTENT(in):: xk
384  REAL(KIND(1d0)), INTENT(in):: xch
385  REAL(KIND(1d0)), INTENT(in):: xbo
386 
387  ! output
388  REAL(KIND(1d0)), INTENT(out) :: xa1
389  REAL(KIND(1d0)), INTENT(out) :: xa2
390  REAL(KIND(1d0)), INTENT(out) :: xa3
391  REAL(KIND(1d0)), INTENT(out) :: ats
392  REAL(KIND(1d0)), INTENT(out) :: mts
393  REAL(KIND(1d0)), INTENT(out) :: gamma
394 
395  ! constant
396  REAL(KIND(1d0)), PARAMETER :: sigma = 5.67e-8 ! Stefan-Boltzman
397  REAL(KIND(1d0)), PARAMETER :: pi = atan(1.0)*4 ! Pi
398  REAL(KIND(1d0)), PARAMETER :: omega = 2*pi/(24*60*60) ! augular velocity of Earth
399 
400  ! local variables:
401  REAL(KIND(1d0)) :: beta ! inverse Bowen ratio
402  REAL(KIND(1d0)) :: f, fl, ft ! energy redistribution factors
403  REAL(KIND(1d0)) :: lambda ! thermal diffusivity
404  REAL(KIND(1d0)) :: delta, m, n ! water flux related variables
405  REAL(KIND(1d0)) :: ms, ns ! m, n related
406  REAL(KIND(1d0)) :: ceta, cphi ! phase related temporary variables
407  REAL(KIND(1d0)) :: eta, phi, xlag ! phase related temporary variables
408  REAL(KIND(1d0)) :: xx1, xx2, xx3, xchws ! temporary use
409  ! LOGICAL :: flagGood = .TRUE. ! quality flag, T for good, F for bad
410 
411  ! give fixed values for test
412  ! properties
413  ! xalb = .2
414  ! xemis = .9
415  ! xcp = 1e6
416  ! xk = 1.2
417  ! xch = 4
418  ! xBo = 2
419  ! forcings
420  ! ASd = 400
421  ! mSd = 100
422  ! ATa = 23
423  ! mTa = 23+C2K
424  ! tau = PI/6
425  ! WS = 1
426  ! AH = 0
427  ! Wf = 0
428 
429  ! PRINT*, '********sfc_typ: ',sfc_typ,' start********'
430  ! initialize flagGood as .TRUE.
431  ! flagGood = .TRUE.
432  ! PRINT*, flagGood
433 
434  ! calculate sfc properties related parameters:
435  xchws = xch*mws
436 
437  ! xch = CRA/RA !update bulk trsf. coeff. with RA (aerodyn. res.)
438  IF (abs(xbo) < 0.001) THEN
439  beta = 1/0.001
440  ELSE
441  beta = 1/xbo
442  END IF
443 
444  ! PRINT*, 'beta:', beta
445  f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
446  ! print*, 'xch',xch,'mTa',mTa,'dLu',4*SIGMA*xemis*mTa**3
447  fl = 4*sigma*xemis*mta**3
448  ! print*, 'fL',fL
449  ft = (1 + beta)*xchws
450  ! print*, 'fT',fT
451  lambda = xk/xcp
452  ! print*, 'lambda',lambda
453  delta = sqrt(.5*(mwf**2 + sqrt(mwf**4 + 16*lambda**2*omega**2)))
454  ! print*, 'delta',delta
455  m = (2*lambda)/(delta + mwf)
456  n = delta/omega
457  ! print*, 'm',m,'n',n
458 
459  ! calculate surface temperature related parameters:
460  mts = (msd*(1 - xalb)/f) + mta
461  ms = 1 + xk/(f*m)
462  ns = xk/(f*n)
463  ! print*, 'ms,ns:', ms,ns
464  xx1 = f*sin(tau)*ata
465  ! print*, 'xx1',xx1
466  xx2 = (1 - xalb)*asd + f*cos(tau)*ata
467  ! print*, 'xx2',xx2
468  gamma = atan(ns/ms) + atan(xx1/xx2)
469  ! print*, 'gamma:', gamma
470  ats = -(sin(tau)*ata)/(ns*cos(gamma) - ms*sin(gamma))
471  ! print*, 'ATs:', ATs
472 
473  ! calculate net radiation related parameters:
474  xx1 = (ns*cos(gamma) + sin(gamma) - ms*sin(gamma))*sin(tau)*ata*fl
475  xx2 = (xalb - 1)*(ns*cos(gamma) - ms*sin(gamma))*asd
476  xx3 = (-ms*cos(tau)*sin(tau) + cos(gamma)*(ns*cos(tau) + sin(tau)))*ata*fl
477  xx2 = xx2 - xx3
478  phi = atan(xx1/xx2)
479  xx3 = (ns*cos(gamma) - ms*sin(gamma))
480  xx1 = (1 + sin(gamma)/xx3)*sin(tau)*ata*fl
481  xx2 = (xalb - 1)*asd - (cos(tau) + cos(gamma)*sin(tau)/xx3)*ata*fl
482  cphi = sqrt(xx1**2 + xx2**2)
483 
484  ! calculate heat storage related parameters:
485  xx1 = m*cos(gamma) - n*sin(gamma)
486  xx2 = m*sin(gamma) + n*cos(gamma)
487  eta = atan(xx1/xx2)
488  ! if ( eta<0 ) print*, 'lambda,delta,m,n,gamma:', lambda,delta,m,n,gamma
489  xx1 = xk**2*(m**2 + n**2)*ats**2
490  xx2 = m**2*n**2
491  ceta = sqrt(xx1/xx2)
492 
493  ! key phase lag:
494  xlag = eta - phi
495 
496  ! calculate the OHM coeffs.:
497  ! a1:
498  xa1 = (ceta*cos(xlag))/cphi
499 
500  ! a2:
501  xa2 = (ceta*sin(xlag))/(omega*cphi)
502  xa2 = xa2/3600 ! convert the unit from s-1 to h-1
503 
504  ! a3:
505  xa3 = -xa1*(ft/f)*(msd*(1 - xalb)) - mah
506 
507  ! print*, 'ceta,cphi', ceta,cphi
508  ! print*, 'tau,eta,phi,xlag in deg:',tau/pi*180,eta/pi*180,phi/pi*180,xlag/pi*180
509 
510  ! PRINT*, '********sfc_typ: ',sfc_typ,' end********'
511 
allocatearray::xbo
real(kind(1d0)) xbo
Definition: suews_ctrl_const.f95:320
Here is the caller graph for this function:

◆ anohm_coef_water_cal()

subroutine anohm_module::anohm_coef_water_cal ( real(kind(1d0)), intent(in)  ASd,
real(kind(1d0)), intent(in)  mSd,
real(kind(1d0)), intent(in)  ATa,
real(kind(1d0)), intent(in)  mTa,
real(kind(1d0)), intent(in)  tau,
real(kind(1d0)), intent(in)  mWS,
real(kind(1d0)), intent(in)  mWF,
real(kind(1d0)), intent(in)  mAH,
real(kind(1d0)), intent(in)  xalb,
real(kind(1d0)), intent(in)  xemis,
real(kind(1d0)), intent(in)  xcp,
real(kind(1d0)), intent(in)  xk,
real(kind(1d0)), intent(in)  xch,
real(kind(1d0)), intent(in)  xBo,
real(kind(1d0)), intent(in)  xeta,
real(kind(1d0)), intent(in)  xmu,
real(kind(1d0)), intent(out)  xa1,
real(kind(1d0)), intent(out)  xa2,
real(kind(1d0)), intent(out)  xa3,
real(kind(1d0)), intent(out)  ATs,
real(kind(1d0)), intent(out)  mTs,
real(kind(1d0)), intent(out)  gamma 
)

a wrapper for retrieving AnOHM coefficients of water body

Returns
  1. OHM coefficients of a given surface type: a1, a2 and a3 : NB: this SUBROUTINE hasn't been well tested.
Parameters
[in]asddaily amplitude of solar radiation [W m-2]
[in]msddaily mean solar radiation [W m-2]
[in]atadaily amplitude of air temperature [K]
[in]mtadaily mean air temperature [K]
[in]tauphase lag between Sd and Ta (Ta-Sd) [rad]
[in]mwsdaily mean wind speed [m s-1]
[in]mwfdaily mean underground moisture flux [m3 s-1 m-2]
[in]mahdaily mean anthropogenic heat flux [W m-2]
[in]xalbalbedo [-]
[in]xemisemissivity [-]
[in]xcpheat capacity [J m-3 K-1]
[in]xkthermal conductivity [W m-1 K-1]
[in]xchbulk transfer coef [J m-3 K-1]
[in]xboBowen ratio [-]
[in]xetaeffective absorption fraction [-]
[in]xmueffective absorption coefficient [m-1]
[out]xa1AnOHM coefficients of grid [-]
[out]xa2AnOHM coefficients of grid [h]
[out]xa3AnOHM coefficients of grid [W m-2]
[out]atsdaily amplitude of surface temperature [K]
[out]mtsdaily mean of surface temperature [K]
[out]gammaphase difference between Ts and Sd [K]

Definition at line 524 of file suews_phys_anohm.f95.

Referenced by anohm_coef(), and anohm_xts().

524 
525  IMPLICIT NONE
526 
527  ! input: forcing scales
528  REAL(KIND(1d0)), INTENT(in):: asd
529  REAL(KIND(1d0)), INTENT(in):: msd
530  REAL(KIND(1d0)), INTENT(in):: ata
531  REAL(KIND(1d0)), INTENT(in):: mta
532  REAL(KIND(1d0)), INTENT(in):: tau
533  REAL(KIND(1d0)), INTENT(in):: mws
534  REAL(KIND(1d0)), INTENT(in):: mwf
535  REAL(KIND(1d0)), INTENT(in):: mah
536 
537  ! input: sfc properties
538  REAL(KIND(1d0)), INTENT(in):: xalb
539  REAL(KIND(1d0)), INTENT(in):: xemis
540  REAL(KIND(1d0)), INTENT(in):: xcp
541  REAL(KIND(1d0)), INTENT(in):: xk
542  REAL(KIND(1d0)), INTENT(in):: xch
543  REAL(KIND(1d0)), INTENT(in):: xbo
544  REAL(KIND(1d0)), INTENT(in):: xeta
545  REAL(KIND(1d0)), INTENT(in):: xmu
546 
547  ! output
548  REAL(KIND(1d0)), INTENT(out) :: xa1
549  REAL(KIND(1d0)), INTENT(out) :: xa2
550  REAL(KIND(1d0)), INTENT(out) :: xa3
551  REAL(KIND(1d0)), INTENT(out) :: ats
552  REAL(KIND(1d0)), INTENT(out) :: mts
553  REAL(KIND(1d0)), INTENT(out) :: gamma
554 
555  ! constant
556  REAL(KIND(1d0)), PARAMETER :: sigma = 5.67e-8 ! Stefan-Boltzman
557  REAL(KIND(1d0)), PARAMETER :: pi = atan(1.0)*4 ! Pi
558  REAL(KIND(1d0)), PARAMETER :: omega = 2*pi/(24*60*60) ! augular velocity of Earth
559 
560  ! local variables:
561  REAL(KIND(1d0)) :: beta ! inverse Bowen ratio
562  REAL(KIND(1d0)) :: f, fl, ft ! energy redistribution factors
563  REAL(KIND(1d0)) :: lambda, calb ! temporary use
564  REAL(KIND(1d0)) :: delta ! sfc. temperature related variables
565  ! REAL(KIND(1d0)) :: delta,m,n ! sfc. temperature related variables
566  REAL(KIND(1d0)) :: xm, xn ! m, n related
567  ! REAL(KIND(1d0)) :: gamma ! phase lag scale
568  REAL(KIND(1d0)) :: phi ! phase lag scale
569  ! REAL(KIND(1d0)) :: ATs,mTs ! surface temperature amplitude
570  REAL(KIND(1d0)) :: czeta, ctheta ! phase related temporary variables
571  REAL(KIND(1d0)) :: zeta, theta, xlag ! phase related temporary variables
572  REAL(KIND(1d0)) :: xx1, xx2, xx3 ! temporary use
573  REAL(KIND(1d0)) :: kappa ! temporary use
574  REAL(KIND(1d0)) :: dtau, dpsi, dphi ! temporary use
575  REAL(KIND(1d0)) :: cdtau, cdpsi, cdphi ! temporary use
576  REAL(KIND(1d0)) :: xxt, xxkappa, xxdltphi, xchws ! temporary use
577  ! LOGICAL :: flagGood = .TRUE. ! quality flag, T for good, F for bad
578 
579  ! ====not used====
580  REAL(KIND(1d0)) :: dummy
581  dummy = mah + mwf
582  ! ====not used====
583 
584  ! calculate sfc properties related parameters:
585  xm = xk*xmu**2
586  xn = xcp*omega
587  phi = atan(xn/xm)
588  kappa = sqrt(xcp*omega/(2*xk))
589 
590  ! mWS=SUM(WS, dim=1, mask=(WS>0))/SIZE(WS, dim=1)
591  xchws = xch*mws
592  beta = 1/xbo
593  f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
594  fl = 4*sigma*xemis*mta**3
595  ft = (1 + beta)*xchws
596 
597  calb = 1 - xalb
598 
599  lambda = sqrt(xm**2 + xn**2)
600 
601  dtau = atan(xk*kappa/(f + xk*kappa))
602  dpsi = atan((xk - xmu)/kappa)
603  dphi = atan(kappa*(f + xk*xmu)/(f*(kappa - xmu) + xk*kappa*(2*kappa - xmu)))
604  cdtau = sqrt((xk*kappa)**2 + (f + xk*kappa)**2)
605  cdpsi = sqrt((xk - xmu)**2 + kappa**2)
606  cdphi = cdtau*cdpsi
607 
608  ! calculate surface temperature related parameters:
609  ! daily mean:
610  mts = (msd*(1 - xalb + xeta)/f) + mta
611  ! amplitude:
612  xx1 = (xk*xeta*xmu*calb*asd*cdpsi)**2
613  xx2 = 2*lambda*sqrt(xx1)*(calb*asd*sin(phi - dpsi) + f*ata*sin(tau + phi - dpsi))
614  xx3 = lambda**2*((calb*asd + cos(tau)*f*ata)**2 + (sin(tau)*f*ata)**2)
615  ats = 1/(cdtau*lambda)*sqrt(xx1 + xx2 + xx3)
616  ! phase lag:
617  xx1 = (xk*kappa*calb*asd + cdtau*f*ata*sin(tau + dtau))*lambda &
618  + (xk*xeta*xmu*calb*asd*cdphi*sin(phi + dphi))
619  xx2 = ((f + xk*kappa)*calb*asd - cdtau*f*ata*cos(tau + dtau))*lambda &
620  - xk*xeta*xmu*calb*asd*cdphi*cos(phi + dphi)
621  delta = atan(xx1/xx2)
622  gamma = delta
623 
624  ! calculate net radiation related parameters:
625  ! phase lag:
626  xx1 = fl*(ats*sin(delta) - ata*sin(tau))
627  xx2 = calb*asd - fl*(ats*cos(delta) - ata*cos(tau))
628  theta = atan(xx1/xx2)
629  ! amplitude:
630  ctheta = sqrt(xx1**2 + xx2**2)
631 
632  ! calculate heat storage related parameters:
633  ! scales:
634  xxt = sqrt(2.)*kappa*lambda*ats
635  xxkappa = cdpsi*xeta*xmu*asd
636  xxdltphi = cos(delta)*sin(dpsi)*cos(phi) - sin(delta)*cos(dpsi)*sin(phi)
637  ! phase lag:
638  xx1 = xxt*sin(pi/4 - delta) + xxkappa*sin(phi + dpsi)
639  xx2 = xxt*sin(pi/4 + delta) - xxkappa*sin(phi - dpsi)
640  zeta = atan(xx1/xx2)
641  ! amplitude:
642  xx1 = 2*sqrt(2.)*xxkappa*xxt*xxdltphi
643  xx2 = (1 - cos(2*dpsi)*cos(2*phi))*xxkappa**2
644  xx3 = xxt**2
645  czeta = xk/lambda*sqrt(xx1 + xx2 + xx3)
646 
647  ! calculate the OHM coeffs.:
648  xlag = zeta - theta
649  ! a1:
650  xa1 = (czeta*cos(xlag))/ctheta
651 
652  ! write(*,*) 'ceta,xlag,cphi:', ceta,xlag,cphi
653  ! a2:
654  xa2 = (czeta*sin(xlag))/(omega*ctheta)
655  xa2 = xa2/3600 ! convert the unit from s-1 to h-1
656 
657  ! a3:
658  xa3 = msd*(xalb - 1)*(xeta + (ft - fl*xeta)/f*xa1)
659 
660  ! print*, 'sfc_typ_water:', sfc_typ
661  ! print*, 'a1,a2,a3:', xa1,xa2,xa3
662 
allocatearray::xbo
real(kind(1d0)) xbo
Definition: suews_ctrl_const.f95:320
Here is the caller graph for this function:

◆ anohm_fc()

subroutine anohm_module::anohm_fc ( integer, intent(in)  xid,
real(kind(1d0)), dimension(:, :), intent(in)  MetForcingData_grid,
integer, intent(in)  EmissionsMethod,
real(kind(1d0)), intent(in)  qf,
real(kind(1d0)), intent(out)  ASd,
real(kind(1d0)), intent(out)  mSd,
real(kind(1d0)), intent(out)  tSd,
real(kind(1d0)), intent(out)  ATa,
real(kind(1d0)), intent(out)  mTa,
real(kind(1d0)), intent(out)  tTa,
real(kind(1d0)), intent(out)  tau,
real(kind(1d0)), intent(out)  mWS,
real(kind(1d0)), intent(out)  mWF,
real(kind(1d0)), intent(out)  mAH 
)
Parameters
[in]qfanthropogenic heat flux [W m-2]
[out]asddaily amplitude of solar radiation [W m-2]
[out]msddaily mean solar radiation [W m-2]
[out]tsdlocal peaking time of solar radiation [hr]
[out]atadaily amplitude of air temperature [degC]
[out]mtadaily mean air temperature [degC]
[out]ttalocal peaking time of air temperature [hour]
[out]tauphase lag between Sd and Ta (Ta-Sd) [rad]
[out]mwsdaily mean wind speed [m s-1]
[out]mwfdaily mean underground moisture flux [m3 s-1 m-2]
[out]mahdaily mean anthropogenic heat flux [W m-2]

Definition at line 670 of file suews_phys_anohm.f95.

References anohm_fccal(), and anohm_fcload().

Referenced by anohm_coef().

670 
671  IMPLICIT NONE
672 
673  ! input
674  INTEGER, INTENT(in):: xid
675  INTEGER, INTENT(in):: emissionsmethod
676  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :) ::metforcingdata_grid
677  REAL(KIND(1d0)), INTENT(in):: qf
678 
679  ! output
680  REAL(KIND(1d0)), INTENT(out):: asd
681  REAL(KIND(1d0)), INTENT(out):: msd
682  REAL(KIND(1d0)), INTENT(out):: tsd
683  REAL(KIND(1d0)), INTENT(out):: ata
684  REAL(KIND(1d0)), INTENT(out):: mta
685  REAL(KIND(1d0)), INTENT(out):: tta
686  REAL(KIND(1d0)), INTENT(out):: tau
687  REAL(KIND(1d0)), INTENT(out):: mws
688  REAL(KIND(1d0)), INTENT(out):: mwf
689  REAL(KIND(1d0)), INTENT(out):: mah
690 
691  ! forcings:
692  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: sd
693  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: ta
694  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: rh
695  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: pres
696  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: ws ! wind speed [m s-1]
697  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: wf ! water flux density [m3 s-1 m-2]
698  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: ah ! anthropogenic heat [W m-2]
699  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: thr ! local time [hr]
700 
701  ! load forcing variables:
702  CALL anohm_fcload(xid, metforcingdata_grid, emissionsmethod, qf, sd, ta, rh, pres, ws, wf, ah, thr)
703  ! calculate forcing scales for AnOHM:
704  CALL anohm_fccal(sd, ta, ws, wf, ah, thr, asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
705 
706  ! CALL r8vec_print(SIZE(sd, dim=1),sd,'Sd')
707  ! PRINT*, ASd,mSd,tSd
708 
709  ! CALL r8vec_print(SIZE(ta, dim=1),ta,'Ta')
710  ! PRINT*, ATa,mTa,tTa
711 
allocatearray::metforcingdata_grid
real(kind(1d0)), dimension(:, :), allocatable metforcingdata_grid
Definition: suews_ctrl_const.f95:116
estm_data::ws
real(kind(1d0)) ws
Definition: suews_phys_estm.f95:85
data_in::emissionsmethod
integer emissionsmethod
Definition: suews_ctrl_const.f95:977
Here is the call graph for this function:
Here is the caller graph for this function:

◆ anohm_fccal()

subroutine anohm_module::anohm_fccal ( real(kind(1d0)), dimension(:), intent(in)  Sd,
real(kind(1d0)), dimension(:), intent(in)  Ta,
real(kind(1d0)), dimension(:), intent(in)  WS,
real(kind(1d0)), dimension(:), intent(in)  WF,
real(kind(1d0)), dimension(:), intent(in)  AH,
real(kind(1d0)), dimension(:), intent(in)  tHr,
real(kind(1d0)), intent(out)  ASd,
real(kind(1d0)), intent(out)  mSd,
real(kind(1d0)), intent(out)  tSd,
real(kind(1d0)), intent(out)  ATa,
real(kind(1d0)), intent(out)  mTa,
real(kind(1d0)), intent(out)  tTa,
real(kind(1d0)), intent(out)  tau,
real(kind(1d0)), intent(out)  mWS,
real(kind(1d0)), intent(out)  mWF,
real(kind(1d0)), intent(out)  mAH 
)

calculate the key parameters of a sinusoidal curve for AnOHM forcings i.e., a, b, c in a*Sin(Pi/12*t+b)+c

Parameters
[in]sdincoming shortwave radiation [W m-2]
[in]taair temperature [degC]
[in]wswind speed [m s-1]
[in]wfwater flux density [m3 s-1 m-2]
[in]ahanthropogenic heat [W m-2]
[in]thrtime [hr]
[out]asddaily amplitude of solar radiation [W m-2]
[out]msddaily mean solar radiation [W m-2]
[out]tsdlocal peaking time of solar radiation [hr]
[out]atadaily amplitude of air temperature [degC]
[out]mtadaily mean air temperature [degC]
[out]ttalocal peaking time of air temperature [hr]
[out]tauphase lag between Sd and Ta (Ta-Sd) [rad]
[out]mwsdaily mean wind speed [m s-1]
[out]mwfdaily mean underground moisture flux [m3 s-1 m-2]
[out]mahdaily mean anthropogenic heat flux [W m-2]

Definition at line 811 of file suews_phys_anohm.f95.

References anohm_shapefit(), and r8vec_print().

Referenced by anohm_fc().

811  IMPLICIT NONE
812 
813  ! input
814  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: sd
815  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: ta
816  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: ws
817  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: wf
818  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: ah
819  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: thr
820 
821  ! output
822  REAL(KIND(1d0)), INTENT(out):: asd
823  REAL(KIND(1d0)), INTENT(out):: msd
824  REAL(KIND(1d0)), INTENT(out):: tsd
825  REAL(KIND(1d0)), INTENT(out):: ata
826  REAL(KIND(1d0)), INTENT(out):: mta
827  REAL(KIND(1d0)), INTENT(out):: tta
828  REAL(KIND(1d0)), INTENT(out):: tau
829  REAL(KIND(1d0)), INTENT(out):: mws
830  REAL(KIND(1d0)), INTENT(out):: mwf
831  REAL(KIND(1d0)), INTENT(out):: mah
832 
833  ! constant
834  REAL(KIND(1d0)), PARAMETER :: pi = atan(1.0)*4 ! Pi
835  REAL(KIND(1d0)), PARAMETER :: c2k = 273.15 ! degC to K
836 
837  ! local variables:
838  REAL(KIND(1d0)), ALLOCATABLE :: thrday(:) ! daytime tHr when Sd>0
839  REAL(KIND(1d0)), ALLOCATABLE :: selx(:) ! daytime sublist of met variable when Sd>0
840 
841  ! REAL(KIND(1d0)) :: xx ! temporary use
842  INTEGER :: err, lenday
843  LOGICAL, DIMENSION(:), ALLOCATABLE::sdmask
844 
845  ALLOCATE (sdmask(SIZE(sd, dim=1)), stat=err)
846  IF (err /= 0) print *, "SdMask: Allocation request denied"
847  sdmask = sd > 5
848  lenday = count(sdmask)
849 
850  ! CALL r8vec_print(24,Sd,'Sd')
851  ! CALL r8vec_print(24,tHr,'tHr')
852  ALLOCATE (thrday(lenday), stat=err)
853  IF (err /= 0) print *, "tHrDay: Allocation request denied"
854  thrday = pack(thr, mask=sdmask)
855 
856  ALLOCATE (selx(lenday), stat=err)
857  IF (err /= 0) print *, "selX: Allocation request denied"
858 
859  ! calculate sinusoidal scales of Sd:
860  ! PRINT*, 'Calc. Sd...'
861  selx = pack(sd, mask=sdmask)
862  asd = (maxval(selx) - minval(selx))/2
863  msd = sum(selx)/lenday
864  tsd = 12
865  CALL anohm_shapefit(thrday, selx, asd, msd, tsd)
866  ! CALL r8vec_print(lenDay,selX,'Sd Day:')
867  ! modify ill-shaped days to go through
868  ! IF ( ASd < 0 .OR. tSd > 15) THEN
869  ! ! ASd = abs(ASd)
870  ! ! tSd = 12 ! assume Sd peaks at 12:00LST
871  ! CALL r8vec_print(lenDay,tHrDay,'tHrDay:')
872  ! CALL r8vec_print(lenDay,selX,'Sd Day:')
873  ! PRINT*, 'ASd:', ASd
874  ! PRINT*, 'mSd:', mSd
875  ! PRINT*, 'tSd:', tSd
876  ! END IF
877  ! PRINT*, 'Sd Day:', selX
878 
879  ! add corrections on mSd to fix the overestimated a3 during cold periods
880  IF (msd + asd < maxval(selx)) THEN
881  ! PRINT*, 'Sd max:', MAXVAL(selX)
882  ! PRINT*, 'ASd:', ASd
883  ! PRINT*, 'mSd before:', mSd
884  msd = maxval(selx) - asd
885  ! PRINT*, 'mSd after:', mSd
886  ! PRINT*,''
887  END IF
888 
889  ! calculate sinusoidal scales of Ta:
890  ! PRINT*, 'Calc. Ta...'
891  selx = pack(ta, mask=sdmask)
892  ata = (maxval(selx) - minval(selx))/2
893  mta = sum(selx)/lenday
894  tta = 12
895  CALL anohm_shapefit(thrday, selx, ata, mta, tta)
896  ! CALL r8vec_print(lenDay,selX,'Ta Day:')
897  IF (mta < 60) mta = mta + c2k ! correct the Celsius to Kelvin
898  ! modify ill-shaped days to go through
899  IF (ata < 0) THEN
900  ! ATa = abs(ATa)
901  ! tTa = 14 ! assume Ta peaks at 14:00LST
902  CALL r8vec_print(lenday, selx, 'Ta Day:')
903  print *, 'ATa:', ata
904  print *, 'mTa:', mta
905  print *, 'tTa:', tta
906  END IF
907  ! PRINT*, 'Ta:', Ta(10:16)
908 
909  ! calculate the phase lag between Sd and Ta:
910  tau = (tta - tsd)/24*2*pi
911  ! PRINT*, 'tau:', tau
912 
913  ! calculate the mean values:
914  selx = pack(ws, mask=sdmask)
915  mws = sum(selx)/lenday ! mean value of WS
916 
917  selx = pack(wf, mask=sdmask)
918  mwf = sum(selx)/lenday ! mean value of WF
919 
920  selx = pack(ah, mask=sdmask)
921  mah = sum(selx)/lenday ! mean value of AH
922  ! PRINT*, 'mWS:', mWS
923  ! print*, 'mWF:', mWF
924  ! print*, 'mAH:', mAH
925  IF (ALLOCATED(sdmask)) DEALLOCATE (sdmask, stat=err)
926  IF (err /= 0) print *, "SdMask: Deallocation request denied"
927 
928  IF (ALLOCATED(thrday)) DEALLOCATE (thrday, stat=err)
929  IF (err /= 0) print *, "tHrDay: Deallocation request denied"
930 
931  IF (ALLOCATED(selx)) DEALLOCATE (selx, stat=err)
932  IF (err /= 0) print *, "selX: Deallocation request denied"
933 
estm_data::ws
real(kind(1d0)) ws
Definition: suews_phys_estm.f95:85
r8vec_print
subroutine r8vec_print(n, a, title)
Definition: suews_util_minpack.f95:5501
Here is the call graph for this function:
Here is the caller graph for this function:

◆ anohm_fcload()

subroutine anohm_module::anohm_fcload ( integer, intent(in)  xid,
real(kind(1d0)), dimension(:, :), intent(in)  MetForcingData_grid,
integer, intent(in)  EmissionsMethod,
real(kind(1d0)), intent(in)  qf,
real(kind(1d0)), dimension(:), intent(out), allocatable  Sd,
real(kind(1d0)), dimension(:), intent(out), allocatable  Ta,
real(kind(1d0)), dimension(:), intent(out), allocatable  RH,
real(kind(1d0)), dimension(:), intent(out), allocatable  pres,
real(kind(1d0)), dimension(:), intent(out), allocatable  WS,
real(kind(1d0)), dimension(:), intent(out), allocatable  WF,
real(kind(1d0)), dimension(:), intent(out), allocatable  AH,
real(kind(1d0)), dimension(:), intent(out), allocatable  tHr 
)

load forcing series for AnOHM_FcCal

Parameters
[in]xidday of year
[in]emissionsmethodAnthropHeat option
[in]metforcingdata_gridmet forcing array of grid
[in]qfanthropogenic heat flux [W m-2]
[out]sdincoming solar radiation [W m-2]
[out]taair temperature [degC]
[out]rhrelative humidity [%]
[out]presatmospheric pressure [mbar]
[out]wswind speed [m s-1]
[out]wfwater flux density [m3 s-1 m-2]
[out]ahanthropogenic heat [W m-2]
[out]thrlocal time [hr]

Definition at line 720 of file suews_phys_anohm.f95.

Referenced by anohm_coef(), and anohm_fc().

720 
721  IMPLICIT NONE
722 
723  ! input
724  INTEGER, INTENT(in):: xid
725  INTEGER, INTENT(in):: emissionsmethod
726  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :) ::metforcingdata_grid
727  REAL(KIND(1d0)), INTENT(in):: qf
728 
729  ! output
730  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: sd
731  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: ta
732  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: rh
733  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: pres
734  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: ws
735  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: wf
736  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: ah
737  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: thr
738 
739  ! local variables:
740  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE :: submet ! subset array of daytime series
741 
742  INTEGER :: err
743  INTEGER :: lenmetdata, nvar
744 
745  LOGICAL, ALLOCATABLE :: metmask(:)
746 
747  ! construct mask
748  IF (ALLOCATED(metmask)) DEALLOCATE (metmask, stat=err)
749  ALLOCATE (metmask(SIZE(metforcingdata_grid, dim=1)))
750  metmask = (metforcingdata_grid(:, 2) == xid & ! day=xid
751  .AND. metforcingdata_grid(:, 4) == 0)! tmin=0
752 
753  ! determine the length of subset
754  lenmetdata = count(metmask)
755 
756  ! construct array for time and met variables
757  nvar = 8! number of variables to retrieve
758  ! print*, 'good 1'
759  ! allocate subMet:
760  IF (ALLOCATED(submet)) DEALLOCATE (submet, stat=err)
761  ALLOCATE (submet(lenmetdata, nvar))
762  submet = reshape(pack(metforcingdata_grid(:, (/3, & !time: hour
763  15, 12, 11, 13, 10, 12, 9/)), &! met: Sd, Ta, RH, pres, WS, WF, AH
764  ! NB: WF not used: a placeholder
765  spread(metmask, dim=2, ncopies=nvar)), & ! replicate mask vector to 2D array
766  (/lenmetdata, nvar/)) ! convert to target shape
767 
768  ! re-allocate arrays as their sizes may change during passing
769  IF (ALLOCATED(thr)) DEALLOCATE (thr, stat=err)
770  ALLOCATE (thr(lenmetdata))
771  IF (ALLOCATED(sd)) DEALLOCATE (sd, stat=err)
772  ALLOCATE (sd(lenmetdata))
773  IF (ALLOCATED(ta)) DEALLOCATE (ta, stat=err)
774  ALLOCATE (ta(lenmetdata))
775  IF (ALLOCATED(rh)) DEALLOCATE (rh, stat=err)
776  ALLOCATE (rh(lenmetdata))
777  IF (ALLOCATED(pres)) DEALLOCATE (pres, stat=err)
778  ALLOCATE (pres(lenmetdata))
779  IF (ALLOCATED(ws)) DEALLOCATE (ws, stat=err)
780  ALLOCATE (ws(lenmetdata))
781  IF (ALLOCATED(wf)) DEALLOCATE (wf, stat=err)
782  ALLOCATE (wf(lenmetdata))
783  IF (ALLOCATED(ah)) DEALLOCATE (ah, stat=err)
784  ALLOCATE (ah(lenmetdata))
785 
786  ! load the sublist into forcing variables
787  thr = submet(:, 1)! time in hour
788  sd = submet(:, 2)
789  ta = submet(:, 3)
790  rh = submet(:, 4)
791  pres = submet(:, 5)
792  ws = submet(:, 6)
793  wf = 0 ! set as 0 for the moment
794  IF (emissionsmethod == 0) THEN
795  ah = submet(:, 8) ! read in from MetForcingData_grid,
796  ELSE
797  ! AH = 0 ! temporarily change to zero;
798  ah = qf ! use modelled value
799  ! AH = mAH_grids(xid-1,xgrid)
800  END IF
801 
allocatearray::metforcingdata_grid
real(kind(1d0)), dimension(:, :), allocatable metforcingdata_grid
Definition: suews_ctrl_const.f95:116
estm_data::ws
real(kind(1d0)) ws
Definition: suews_phys_estm.f95:85
data_in::emissionsmethod
integer emissionsmethod
Definition: suews_ctrl_const.f95:977
Here is the caller graph for this function:

◆ anohm_shapefit()

subroutine anohm_module::anohm_shapefit ( real(kind(1d0)), dimension(:), intent(in)  tHr,
real(kind(1d0)), dimension(:), intent(in)  obs,
real(kind(1d0)), intent(out)  amp,
real(kind(1d0)), intent(out)  mean,
real(kind(1d0)), intent(out)  tpeak 
)

calculate the key parameters of a sinusoidal curve for AnOHM forcings i.e., a, b, c in a*Sin(Pi/12*t+b)+c, where t is in hour

Parameters
[in]thrtime in hour
[in]obsobservation
[out]ampamplitude
[out]meanaverage
[out]tpeakpeaking time (h)

Definition at line 943 of file suews_phys_anohm.f95.

References fsin(), and lmdif1().

Referenced by anohm_fccal().

943  IMPLICIT NONE
944 
945  ! input
946  REAL(KIND(1d0)), INTENT(in) :: thr(:)
947  REAL(KIND(1d0)), INTENT(in) :: obs(:)
948 
949  ! output
950  REAL(KIND(1d0)), INTENT(out) :: amp
951  REAL(KIND(1d0)), INTENT(out) :: mean
952  REAL(KIND(1d0)), INTENT(out) :: tpeak
953 
954  INTEGER :: m, n, info, err ! temporary use
955 
956  ! EXTERNAL fSin
957  REAL(KIND(1d0)), ALLOCATABLE:: fvec(:), x(:)
958 
959  REAL(KIND(1d0)):: tol = 0.00001d+00 ! tolerance
960 
961  n = 3 ! number of parameters to determine
962  m = SIZE(thr, dim=1) ! number of observation pairs
963  ! PRINT*, 'm',m
964 
965  ALLOCATE (fvec(m), stat=err)
966  IF (err /= 0) print *, "fvec: Allocation request denied"
967 
968  ALLOCATE (x(n), stat=err)
969  IF (err /= 0) print *, "x: Allocation request denied"
970 
971  ! initial guess for fitting
972  x = (/mean, amp, tpeak/)
973 
974  ! use minpack subroutine 'lmstr1' to fit the sinusoidal form
975  CALL lmdif1(fsin, m, n, x, thr, obs, fvec, tol, info)
976 
977  ! x = (/mean,amp,delta/), see subroutine 'fSin' for the meaning of this vector
978  mean = x(1)
979  amp = x(2)
980  tpeak = x(3) + 6 ! when t = delta + 6 (hr of LST), t is tpeak.
981 
982  ! adjust fitted parameters to physical range:
983  ! amp>0
984  IF (amp < 0) THEN
985  ! PRINT*, ''
986  ! PRINT*, 'before:'
987  ! PRINT*, amp,tpeak,mean
988  ! CALL r8vec_print(SIZE(tHR, dim=1),tHR,'tHR')
989  ! CALL r8vec_print(SIZE(obs, dim=1),obs,'obs')
990  amp = abs(amp)
991  tpeak = x(3) - 12 + 6 + 24
992  tpeak = mod(tpeak, 24.)
993  ! PRINT*, 'after:'
994  ! PRINT*, amp,tpeak,mean
995  END IF
996  tpeak = mod(tpeak, 24.)
997 
lmdif1
subroutine lmdif1(fcn, m, n, x, xdat, ydat, fvec, tol, info)
Definition: suews_util_minpack.f95:3478
Here is the call graph for this function:
Here is the caller graph for this function:

◆ anohm_xts()

subroutine anohm_module::anohm_xts ( integer, intent(in)  sfc_typ,
real(kind(1d0)), intent(in)  ASd,
real(kind(1d0)), intent(in)  mSd,
real(kind(1d0)), intent(in)  ATa,
real(kind(1d0)), intent(in)  mTa,
real(kind(1d0)), intent(in)  tau,
real(kind(1d0)), intent(in)  mWS,
real(kind(1d0)), intent(in)  mWF,
real(kind(1d0)), intent(in)  mAH,
real(kind(1d0)), intent(in)  xalb,
real(kind(1d0)), intent(in)  xemis,
real(kind(1d0)), intent(in)  xcp,
real(kind(1d0)), intent(in)  xk,
real(kind(1d0)), intent(in)  xch,
real(kind(1d0)), intent(in)  xBo,
real(kind(1d0)), intent(in)  tSd,
real(kind(1d0)), intent(in)  xTHr,
real(kind(1d0)), intent(out)  xTs 
)

calculate the surface temperature related parameters (ATs, mTs, gamma) based on forcings and sfc.

conditions

Returns
xTs surface temperature at local time
Parameters
[in]sfc_typsurface type (land: 1-6, water: 7)
[in]asddaily amplitude of solar radiation [W m-2]
[in]msddaily mean solar radiation [W m-2]
[in]atadaily amplitude of air temperature [K]
[in]mtadaily mean air temperature [K]
[in]tauphase lag between Sd and Ta (Ta-Sd) [rad]
[in]mwsdaily mean wind speed [m s-1]
[in]mwfdaily mean underground moisture flux [m3 s-1 m-2]
[in]mahdaily mean anthropogenic heat flux [W m-2]
[in]xalbalbedo [-]
[in]xemisemissivity [-]
[in]xcpheat capacity [J m-3 K-1]
[in]xkthermal conductivity [W m-1 K-1]
[in]xchbulk transfer coef [J m-3 K-1]
[in]xboBowen ratio [-]
[in]tsdlocal peaking time of Sd, hour
[in]xthrlocal time to calculate Ts, hour
[out]xtssurface temperature at xTHr(hr)

Definition at line 297 of file suews_phys_anohm.f95.

References anohm_coef_land_cal(), and anohm_coef_water_cal().

Referenced by fcnbo().

297 
298  IMPLICIT NONE
299  ! input:
300  INTEGER, INTENT(in):: sfc_typ
301 
302  ! input: forcing scales
303  REAL(KIND(1d0)), INTENT(in):: asd
304  REAL(KIND(1d0)), INTENT(in):: msd
305  REAL(KIND(1d0)), INTENT(in):: ata
306  REAL(KIND(1d0)), INTENT(in):: mta
307  REAL(KIND(1d0)), INTENT(in):: tau
308  REAL(KIND(1d0)), INTENT(in):: mws
309  REAL(KIND(1d0)), INTENT(in):: mwf
310  REAL(KIND(1d0)), INTENT(in):: mah
311 
312  ! input: sfc properties
313  REAL(KIND(1d0)), INTENT(in):: xalb
314  REAL(KIND(1d0)), INTENT(in):: xemis
315  REAL(KIND(1d0)), INTENT(in):: xcp
316  REAL(KIND(1d0)), INTENT(in):: xk
317  REAL(KIND(1d0)), INTENT(in):: xch
318  REAL(KIND(1d0)), INTENT(in):: xbo
319  REAL(KIND(1d0)):: xeta
320  REAL(KIND(1d0)):: xmu
321 
322  ! input: temporal-related
323  REAL(KIND(1d0)), INTENT(in):: tsd
324  REAL(KIND(1d0)), INTENT(in):: xthr
325 
326  ! output:
327  REAL(KIND(1d0)), INTENT(out) :: xts
328 
329  ! local
330  REAL(KIND(1d0)) :: &
331  xa1, xa2, xa3, &!coefficients
332  ats, mts, gamma !surface temperature related scales by AnOHM
333 
334  ! constant:
335  REAL(KIND(1d0)), PARAMETER :: pi = atan(1.0)*4 ! Pi
336  REAL(KIND(1d0)), PARAMETER :: omega = 2*pi/(24*60*60) ! augular velocity of Earth
337 
338  SELECT CASE (sfc_typ)
339  CASE (1:6) ! land surfaces
340  CALL anohm_coef_land_cal( &
341  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
342  xalb, xemis, xcp, xk, xch, xbo, & ! input: sfc properties
343  xa1, xa2, xa3, ats, mts, gamma) ! output: surface temperature related scales by AnOHM
344 
345  CASE (7) ! water surface
346  ! ! NB:give fixed values for the moment
347  xeta = 0.3
348  xmu = 0.2
349  CALL anohm_coef_water_cal( &
350  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
351  xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, & ! input: sfc properties
352  xa1, xa2, xa3, ats, mts, gamma) ! output
353 
354  END SELECT
355 
356  ! for a local time xTHr (in hour):
357  xts = ats*sin(omega*(xthr - tsd + 6)*3600 - gamma) + mts
358 
allocatearray::xbo
real(kind(1d0)) xbo
Definition: suews_ctrl_const.f95:320
Here is the call graph for this function:
Here is the caller graph for this function:

◆ esat_fn()

real(kind(1d0)) function anohm_module::esat_fn ( real(kind(1d0))  Ta)

calculate saturation vapor pressure (es) at air temperature (Ta) (MRR, 1987)

Parameters
taair temperature [degC]
Returns
saturation vapor pressure [hPa]

Definition at line 1401 of file suews_phys_anohm.f95.

Referenced by qa_fn(), and qsat_fn().

1401  REAL(KIND(1D0))::ta
1402  REAL(KIND(1D0))::esat
1403 
1404  REAL(KIND(1D0)), PARAMETER::a = 6.106
1405  REAL(KIND(1D0)), PARAMETER::b = 17.27
1406  REAL(KIND(1D0)), PARAMETER::c = 237.3
1407 
1408  esat = a*exp(b*ta/(c + ta))
Here is the caller graph for this function:

◆ fcnbo()

subroutine anohm_module::fcnbo ( integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  fvec,
integer(kind=4)  iflag,
integer(kind=4)  m,
real(kind=8), dimension(m)  prms 
)

this fucntion will construct an equaiton for Bo calculation

Definition at line 1171 of file suews_phys_anohm.f95.

References anohm_xts(), qa_fn(), qsat_fn(), and r8vec_print().

Referenced by anohm_bo_cal().

1171 
1172  IMPLICIT NONE
1173  ! Input, external FCN, the name of the user-supplied subroutine which
1174  ! calculates the functions. The routine should have the form:
1175  ! subroutine fcn ( n, x, fvec, iflag )
1176  INTEGER(kind=4):: n
1177  INTEGER(kind=4):: m
1178  INTEGER(kind=4):: iflag
1179  REAL(kind=8):: fvec(n)
1180  REAL(kind=8):: x(n) ! x(1) as unknown Bo
1181  REAL(kind=8):: prms(m) ! prms(i) used for passing parameters
1182 
1183  ! the unknow: Bowen ratio
1184  REAL(kind=8):: xbo
1185 
1186  ! forcing scales
1187  REAL(kind=8):: asd
1188  REAL(kind=8):: msd
1189  REAL(kind=8):: ata
1190  REAL(kind=8):: mta
1191  REAL(kind=8):: tau
1192  REAL(kind=8):: mws
1193  REAL(kind=8):: mwf
1194  REAL(kind=8):: mah
1195 
1196  ! sfc. property scales
1197  REAL(kind=8):: xalb
1198  REAL(kind=8):: xemis
1199  REAL(kind=8):: xcp
1200  REAL(kind=8):: xk
1201  REAL(kind=8):: xch
1202 
1203  ! peaking time of Sd in hour
1204  REAL(kind=8)::tsd
1205 
1206  ! surface moisture status [-]
1207  REAL(kind=8)::xsm
1208 
1209  ! surface type
1210  INTEGER :: sfc_typ
1211 
1212  ! array of daytime series
1213  REAL(kind=8), DIMENSION(:, :), ALLOCATABLE :: dayarray
1214 
1215  ! daytime series
1216  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: sd
1217  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: ta
1218  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: rh
1219  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: pres
1220  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: thr
1221  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: ts ! surface temperature, degC
1222  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: qa ! specific humidity, kg kg-1
1223  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: qs ! Saturated specific humidity at surface,kg kg-1
1224 
1225  ! conversion constant:
1226  REAL(kind=8), PARAMETER:: c2k = 273.15 ! degC to K
1227 
1228  ! thermodynamic values at standard temperature (293.15 K) and pressure (101325 pascals):
1229  REAL(kind=8), PARAMETER:: cp_air = 1006.38 ! specific heat capacity of air, J kg-1 K-1
1230  REAL(kind=8), PARAMETER:: lv_air = 2264.705e3 ! latent heat of vaporization of water, J kg-1
1231 
1232  INTEGER :: lenday, i, err, nvar, nprm
1233 
1234  IF (iflag == 0) THEN
1235  ! only print out the x vector
1236 
1237  WRITE (*, '(a)') ''
1238  CALL r8vec_print(n, x, 'x in fcnBo')
1239 
1240  ELSE IF (iflag == 1) THEN
1241  ! calculate fvec at x:
1242  ! pass unknown:
1243  xbo = x(1)
1244 
1245  ! pass forcing scales:
1246  asd = prms(1)
1247  msd = prms(2)
1248  ata = prms(3)
1249  mta = prms(4)
1250  tau = prms(5)
1251  mws = prms(6)
1252  mwf = prms(7)
1253  mah = prms(8)
1254 
1255  ! pass sfc. property scales:
1256  xalb = prms(9)
1257  xemis = prms(10)
1258  xcp = prms(11)
1259  xk = prms(12)
1260  xch = prms(13)
1261  xsm = min(prms(14), 1.0)
1262 
1263  ! pass tSd:
1264  tsd = prms(15)
1265 
1266  ! pass tSd:
1267  sfc_typ = int(prms(16))
1268 
1269  ! set number of parameters according to the above code
1270  nprm = 16
1271 
1272  ! number of met variables to pass:
1273  nvar = 5
1274  ! length of daytime series
1275  lenday = (m - nprm)/nvar
1276  ! allocate daytime series
1277  ALLOCATE (dayarray(nvar, lenday), stat=err)
1278  IF (err /= 0) print *, "dayArray: Allocation request denied"
1279  dayarray = reshape(prms(nprm + 1:SIZE(prms)), shape=(/nvar, lenday/), order=(/2, 1/))
1280 
1281  ! pass daytime series
1282  ! Sd:
1283  ALLOCATE (sd(lenday), stat=err)
1284  IF (err /= 0) print *, "Sd: Allocation request denied"
1285  sd(:) = dayarray(1, :)
1286  ! Ta:
1287  ALLOCATE (ta(lenday), stat=err)
1288  IF (err /= 0) print *, "Ta: Allocation request denied"
1289  ta(:) = dayarray(2, :)
1290  ! RH:
1291  ALLOCATE (rh(lenday), stat=err)
1292  IF (err /= 0) print *, "RH: Allocation request denied"
1293  rh(:) = dayarray(3, :)
1294  ! pres:
1295  ALLOCATE (pres(lenday), stat=err)
1296  IF (err /= 0) print *, "pres: Allocation request denied"
1297  pres(:) = dayarray(4, :)
1298  ! tHr:
1299  ALLOCATE (thr(lenday), stat=err)
1300  IF (err /= 0) print *, "tHr: Allocation request denied"
1301  thr(:) = dayarray(5, :)
1302 
1303  ! PRINT*, 'n',n
1304  ! PRINT*, 'lenDay',lenDay
1305  ! PRINT*, 'nVar',nVar
1306  ! PRINT*, 'shape of dayArray',SHAPE(dayArray)
1307  ! PRINT*, 'shape of met variables',SHAPE(Ta),SHAPE(RH),SHAPE(pres),SHAPE(tHr)
1308 
1309  ALLOCATE (ts(lenday), stat=err)
1310  IF (err /= 0) print *, "Ts: Allocation request denied"
1311  ALLOCATE (qs(lenday), stat=err)
1312  IF (err /= 0) print *, "qs: Allocation request denied"
1313  ALLOCATE (qa(lenday), stat=err)
1314  IF (err /= 0) print *, "qa: Allocation request denied"
1315 
1316  ! calculate sums of QH and QE series in the daytime
1317  IF (xsm == 0) THEN
1318  xbo = 1000 ! extremely arid
1319  ELSE
1320  ! PRINT*, 'lenDay',lenDay
1321  DO i = 1, lenday, 1
1322  ! calculate surface temperature
1323  CALL anohm_xts( &
1324  sfc_typ, &
1325  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
1326  xalb, xemis, xcp, xk, xch, xbo, & ! input: sfc properties
1327  tsd, & !input: peaking time of Sd in hour
1328  thr(i), &! input: time in hour
1329  ts(i))! output: surface temperature, K
1330 
1331  ! convert K to degC
1332  ts(i) = min(ts(i) - c2k, -40.)
1333 
1334  ! calculate saturation specific humidity
1335  qs(i) = qsat_fn(ts(i), pres(i))
1336 
1337  ! calculate specific humidity
1338  qa(i) = qa_fn(ta(i), rh(i), pres(i))
1339 
1340  ! PRINT*,''
1341  ! PRINT*, 'tHr',tHr(i)
1342  ! PRINT*, 'Sd',Sd(i)
1343  ! PRINT*, 'Ts',Ts(i)
1344  ! PRINT*, 'pres',pres(i)
1345  ! PRINT*, 'qs',qs(i)
1346  ! PRINT*, 'Ta',Ta(i)
1347  ! PRINT*, 'RH',RH(i)
1348  ! PRINT*, 'pres',pres(i)
1349  ! PRINT*, 'qa',qa(i)
1350 
1351  END DO
1352 
1353  ! below for testing:
1354  ! rho_air=1.293, air density, kg m-3
1355  ! RA=60, nominal aerodynamic resistence, s m-1
1356  ! PRINT*,''
1357  ! PRINT*, 'QH:',SUM(cp_air*(Ts-Ta)*1.293/60)/lenDay
1358  ! PRINT*, 'xSM:',xSM
1359  ! PRINT*, 'QE:',SUM(xSM*Lv_air*(qs-qa)*1.293/60)/lenDay
1360  ! PRINT*,''
1361 
1362  xbo = sum(cp_air*(ts - ta))/ & ! sum(QH)
1363  sum(xsm*lv_air*(qs - qa))! sum(QE)
1364 
1365  END IF
1366 
1367  ! f(Bo)-Bo==0
1368  fvec(1) = x(1) - xbo
1369 
1370  ! deallocate arrays
1371  IF (ALLOCATED(dayarray)) DEALLOCATE (dayarray, stat=err)
1372  IF (err /= 0) print *, "dayArray: Deallocation request denied"
1373  IF (ALLOCATED(sd)) DEALLOCATE (sd, stat=err)
1374  IF (err /= 0) print *, "Sd: Deallocation request denied"
1375  IF (ALLOCATED(ta)) DEALLOCATE (ta, stat=err)
1376  IF (err /= 0) print *, "Ta: Deallocation request denied"
1377  IF (ALLOCATED(rh)) DEALLOCATE (rh, stat=err)
1378  IF (err /= 0) print *, "RH: Deallocation request denied"
1379  IF (ALLOCATED(pres)) DEALLOCATE (pres, stat=err)
1380  IF (err /= 0) print *, "pres: Deallocation request denied"
1381  IF (ALLOCATED(thr)) DEALLOCATE (thr, stat=err)
1382  IF (err /= 0) print *, "tHr: Deallocation request denied"
1383  IF (ALLOCATED(ts)) DEALLOCATE (ts, stat=err)
1384  IF (err /= 0) print *, "Ts: Deallocation request denied"
1385  IF (ALLOCATED(qa)) DEALLOCATE (qa, stat=err)
1386  IF (err /= 0) print *, "qa: Deallocation request denied"
1387  IF (ALLOCATED(qs)) DEALLOCATE (qs, stat=err)
1388  IF (err /= 0) print *, "qs: Deallocation request denied"
1389 
1390  END IF
1391  RETURN
1392 
r8vec_print
subroutine r8vec_print(n, a, title)
Definition: suews_util_minpack.f95:5501
allocatearray::xbo
real(kind(1d0)) xbo
Definition: suews_ctrl_const.f95:320
Here is the call graph for this function:
Here is the caller graph for this function:

◆ fsin()

subroutine anohm_module::fsin ( integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  xdat,
real(kind=8), dimension(m)  ydat,
real(kind=8), dimension(m)  fvec,
integer(kind=4)  iflag 
)

sinusoidal function f(t) for fitting: f(t) = mean+amp*Sin(Pi/12(t-delta)) x = (/mean,amp,delta/) contains the fitting parameters

Definition at line 1006 of file suews_phys_anohm.f95.

Referenced by anohm_shapefit().

1006 
1007  IMPLICIT NONE
1008  INTEGER(kind=4) m
1009  INTEGER(kind=4) n
1010 
1011  ! REAL ( kind = 8 ) fjrow(n)
1012  REAL(kind=8) fvec(m), & ! residual vector
1013  xdat(m), ydat(m) ! (x,y) observations for fitting
1014 
1015  INTEGER(kind=4) iflag, i
1016  REAL(kind=8) x(n)
1017  REAL(KIND(1d0)), PARAMETER :: pi = atan(1.0)*4 ! Pi
1018 
1019  IF (iflag == 0) THEN
1020 
1021  WRITE (*, '(a)') ''
1022  DO i = 1, n
1023  WRITE (*, '(g14.6)') x(i)
1024  END DO
1025 
1026  ELSE IF (iflag == 1) THEN
1027  fvec(1:m) = x(1) + x(2)*sin(2*pi/24*(xdat(1:m) - x(3))) - ydat(1:m)
1028 
1029  END IF
1030  RETURN
1031 
Here is the caller graph for this function:

◆ qa_fn()

real(kind(1d0)) function anohm_module::qa_fn ( real(kind(1d0))  Ta,
real(kind(1d0))  RH,
real(kind(1d0))  pres 
)

convert relative humidity (RH) to specific humidity (qa) at air temperature (Ta) and atmospheric pressure (pres)

Parameters
taair temperature [degC]
rhrelative humidity [%]
presatmospheric pressure [hPa]
Returns
saturation specific humidity [kg kg-1]

Definition at line 1434 of file suews_phys_anohm.f95.

References esat_fn().

Referenced by fcnbo().

1434 
1435  REAL(KIND(1D0))::ta
1436  REAL(KIND(1D0))::rh
1437  REAL(KIND(1D0))::ea
1438  REAL(KIND(1D0))::es
1439  REAL(KIND(1D0))::pres
1440  REAL(KIND(1D0))::qa
1441 
1442  REAL(KIND(1D0)), PARAMETER::molar = 0.028965
1443  REAL(KIND(1D0)), PARAMETER::molar_wat_vap = 0.0180153
1444 
1445  es = esat_fn(ta)
1446  ea = es*rh/100
1447  qa = (molar_wat_vap/molar)*ea/pres!(rmh2o/rmair)*ES/PMB
Here is the call graph for this function:
Here is the caller graph for this function:

◆ qsat_fn()

real(kind(1d0)) function anohm_module::qsat_fn ( real(kind(1d0))  Ta,
real(kind(1d0))  pres 
)

calculate saturation specific humidity (qsat) at air temperature (Ta) and atmospheric pressure (pres) (MRR, 1987)

Parameters
taair temperature [degC]
Returns
saturation specific humidity [kg kg-1]
Parameters
presatmospheric pressure [hPa]

Definition at line 1417 of file suews_phys_anohm.f95.

References esat_fn().

Referenced by fcnbo().

1417  REAL(KIND(1D0))::ta
1418  REAL(KIND(1D0))::es
1419  REAL(KIND(1D0))::qsat
1420  REAL(KIND(1D0))::pres
1421 
1422  REAL(KIND(1D0)), PARAMETER::molar = 0.028965
1423  REAL(KIND(1D0)), PARAMETER::molar_wat_vap = 0.0180153
1424 
1425  es = esat_fn(ta)
1426  qsat = (molar_wat_vap/molar)*es/pres!(rmh2o/rmair)*ES/PMB
Here is the call graph for this function:
Here is the caller graph for this function:
Generated on Fri May 15 2020 09:18:05 for SUEWS API Site by   doxygen 1.8.13