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.
 
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.
 
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.
 
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
 
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
 
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
 
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
 
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
 
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
 
subroutine fcnbo (n, x, fvec, iflag, m, prms)
 this fucntion will construct an equaiton for Bo calculation
 
real(kind(1d0)) function esat_fn (ta)
 calculate saturation vapor pressure (es) at air temperature (Ta) (MRR, 1987)
 
real(kind(1d0)) function qsat_fn (ta, pres)
 calculate saturation specific humidity (qsat) at air temperature (Ta) and atmospheric pressure (pres) (MRR, 1987)
 
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)
 

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
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:338
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:262

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

Here is the call 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
hybrd1
subroutine hybrd1(fcn, n, x, fvec, tol, info, m, prms)
Definition suews_util_minpack.f95:1414

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
r8vec_print
subroutine r8vec_print(n, a, title)
Definition suews_util_minpack.f95:5501

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
lmdif1
subroutine lmdif1(fcn, m, n, x, xdat, ydat, fvec, tol, info)
Definition suews_util_minpack.f95:3478

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:
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0