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, qn_av_prev, dqndt_prev, qf, MetForcingData_grid, moist_surf, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, sfr_surf, nsurf, EmissionsMethod, id, Gridiv, qn_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)  qn_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_surf,
integer, intent(in)  nsurf,
integer, intent(in)  EmissionsMethod,
integer, intent(in)  id,
integer, intent(in)  Gridiv,
real(kind(1d0)), intent(out)  qn_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]sfr_surfsurface 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 28 of file suews_phys_anohm.f95.

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_surf(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) :: qn_av_prev
61 REAL(KIND(1D0)), INTENT(out) :: qn_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
84 ! initialize output variables
85 xa1 = 0.1
86 xa2 = 0.2
87 xa3 = 10
88 qs = -999
89 deltaqi = 0 ! NB: as snow part is not implemented within AnOHM yet, this is just a placeholder
90
91 ! to test if the current met block contains enough data for AnOHM
92 ! TODO: more robust selection should be implemented
93 ! daylight hours >= 6
94 idq = count(metforcingdata_grid(:, 2) == id .AND. & ! day of year
95 metforcingdata_grid(:, 4) == 0 .AND. & ! minutes
96 metforcingdata_grid(:, 15) > 0) & ! Sd
97 >= 6
98
99 ! PRINT*, idQ
100 IF (idq) THEN
101 ! given enough data, calculate coefficients of day `id`
102 xid = id
103 id_save = id ! store index of the valid day with enough data
104 ELSE
105 ! otherwise calculate coefficients of yesterday: `id-1`
106 xid = id_save
107 END IF
108
109 DO is = 1, nsurf
110 ! call AnOHM to calculate the coefs.
111 CALL anohm_coef(is, xid, gridiv, metforcingdata_grid, moist_surf, emissionsmethod, qf, & !input
112 alb, emis, cpanohm, kkanohm, chanohm, & ! input
113 xa1(is), xa2(is), xa3(is)) ! output
114 ! print*, 'AnOHM_coef are: ',xa1,xa2,xa3
115 END DO
116
117 ! calculate the areally-weighted OHM coefficients
118 a1 = dot_product(xa1, sfr_surf)
119 a2 = dot_product(xa2, sfr_surf)
120 a3 = dot_product(xa3, sfr_surf)
121
122 ! Calculate radiation part ------------------------------------------------------------
123 qs = -999 !qs = Net storage heat flux [W m-2]
124 IF (qn1 > -999) THEN !qn1 = Net all-wave radiation [W m-2]
125
126 ! Store instantaneous qn1 values for previous hour (qn1_store) and average (qn1_av)
127 ! CALL OHM_dqndt_cal(nsh,qn1,qn1_store,qn1_av_store,dqndt)
128 CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn_av_prev, qn1, dqndt_prev, &
129 qn_av_next, dqndt_next)
130
131 ! Calculate net storage heat flux
132 CALL ohm_qs_cal(qn1, dqndt_prev, a1, a2, a3, qs)
133
134 ELSE
135 CALL errorhint(21, 'SUEWS_AnOHM.f95: bad value for qn found during qs calculation.', qn1, notused, notusedi)
136 END IF
137
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine ohm_qs_cal(qn1, dqndt, a1, a2, a3, qs)
subroutine ohm_dqndt_cal_x(dt, dt_since_start, qn1_av_prev, qn1, dqndt_prev, qn1_av_next, dqndt_next)

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

Referenced by suews_driver::suews_cal_qs().

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 1035 of file suews_phys_anohm.f95.

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

References fcnbo(), and hybrd1().

Referenced by anohm_coef().

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 147 of file suews_phys_anohm.f95.

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

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

Referenced by anohm().

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 361 of file suews_phys_anohm.f95.

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

Referenced by anohm_coef(), and anohm_xts().

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 @Caution: 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 518 of file suews_phys_anohm.f95.

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

Referenced by anohm_coef(), and anohm_xts().

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 665 of file suews_phys_anohm.f95.

668
669 IMPLICIT NONE
670
671 ! input
672 INTEGER, INTENT(in) :: xid
673 INTEGER, INTENT(in) :: EmissionsMethod
674 REAL(KIND(1D0)), INTENT(in), DIMENSION(:, :) :: MetForcingData_grid
675 REAL(KIND(1D0)), INTENT(in) :: qf
676
677 ! output
678 REAL(KIND(1D0)), INTENT(out) :: ASd
679 REAL(KIND(1D0)), INTENT(out) :: mSd
680 REAL(KIND(1D0)), INTENT(out) :: tSd
681 REAL(KIND(1D0)), INTENT(out) :: ATa
682 REAL(KIND(1D0)), INTENT(out) :: mTa
683 REAL(KIND(1D0)), INTENT(out) :: tTa
684 REAL(KIND(1D0)), INTENT(out) :: tau
685 REAL(KIND(1D0)), INTENT(out) :: mWS
686 REAL(KIND(1D0)), INTENT(out) :: mWF
687 REAL(KIND(1D0)), INTENT(out) :: mAH
688
689 ! forcings:
690 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: Sd
691 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: Ta
692 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: RH
693 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: pres
694 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: WS ! wind speed [m s-1]
695 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: WF ! water flux density [m3 s-1 m-2]
696 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: AH ! anthropogenic heat [W m-2]
697 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tHr ! local time [hr]
698
699 ! load forcing variables:
700 CALL anohm_fcload(xid, metforcingdata_grid, emissionsmethod, qf, sd, ta, rh, pres, ws, wf, ah, thr)
701 ! calculate forcing scales for AnOHM:
702 CALL anohm_fccal(sd, ta, ws, wf, ah, thr, asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
703
704 ! CALL r8vec_print(SIZE(sd, dim=1),sd,'Sd')
705 ! PRINT*, ASd,mSd,tSd
706
707 ! CALL r8vec_print(SIZE(ta, dim=1),ta,'Ta')
708 ! PRINT*, ATa,mTa,tTa
709

References anohm_fccal(), and anohm_fcload().

Referenced by anohm_coef().

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 806 of file suews_phys_anohm.f95.

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

References anohm_shapefit(), and r8vec_print().

Referenced by anohm_fc().

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 715 of file suews_phys_anohm.f95.

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

Referenced by anohm_coef(), and anohm_fc().

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 938 of file suews_phys_anohm.f95.

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

References fsin(), and lmdif1().

Referenced by anohm_fccal().

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 288 of file suews_phys_anohm.f95.

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

References anohm_coef_land_cal(), and anohm_coef_water_cal().

Referenced by fcnbo().

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 1398 of file suews_phys_anohm.f95.

1399 REAL(KIND(1D0)) :: Ta
1400 REAL(KIND(1D0)) :: esat
1401
1402 REAL(KIND(1D0)), PARAMETER :: A = 6.106
1403 REAL(KIND(1D0)), PARAMETER :: B = 17.27
1404 REAL(KIND(1D0)), PARAMETER :: C = 237.3
1405
1406 esat = a*exp(b*ta/(c + ta))

Referenced by qa_fn(), and qsat_fn().

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 1168 of file suews_phys_anohm.f95.

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

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

Referenced by anohm_bo_cal().

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 1003 of file suews_phys_anohm.f95.

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

Referenced by anohm_shapefit().

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 1431 of file suews_phys_anohm.f95.

1432
1433 REAL(KIND(1D0)) :: Ta
1434 REAL(KIND(1D0)) :: RH
1435 REAL(KIND(1D0)) :: ea
1436 REAL(KIND(1D0)) :: es
1437 REAL(KIND(1D0)) :: pres
1438 REAL(KIND(1D0)) :: qa
1439
1440 REAL(KIND(1D0)), PARAMETER :: molar = 0.028965
1441 REAL(KIND(1D0)), PARAMETER :: molar_wat_vap = 0.0180153
1442
1443 es = esat_fn(ta)
1444 ea = es*rh/100
1445 qa = (molar_wat_vap/molar)*ea/pres !(rmh2o/rmair)*ES/PMB

References esat_fn().

Referenced by fcnbo().

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 1414 of file suews_phys_anohm.f95.

1415 REAL(KIND(1D0)) :: Ta
1416 REAL(KIND(1D0)) :: es
1417 REAL(KIND(1D0)) :: qsat
1418 REAL(KIND(1D0)) :: pres
1419
1420 REAL(KIND(1D0)), PARAMETER :: molar = 0.028965
1421 REAL(KIND(1D0)), PARAMETER :: molar_wat_vap = 0.0180153
1422
1423 es = esat_fn(ta)
1424 qsat = (molar_wat_vap/molar)*es/pres !(rmh2o/rmair)*ES/PMB

References esat_fn().

Referenced by fcnbo().

Here is the call graph for this function:
Here is the caller graph for this function: