SUEWS API Site
Documentation of SUEWS source code
suews_phys_anohm.f95
Go to the documentation of this file.
1 !========================================================================================
9 ! history:
10 ! 20160301: initial version
11 ! 20170109: updated dqndt calculation in accordance with SUEWS_OHM.f95 (HCW)
12 ! 20170810: revamped structure
13 ! 20170825: improved Bowen calculation
14 !========================================================================================
16 
17  IMPLICIT NONE
18 CONTAINS
19 
20  !========================================================================================
28  SUBROUTINE anohm( &
29  tstep, dt_since_start, &
30  qn1, qn1_av_prev, dqndt_prev, qf, &
31  MetForcingData_grid, moist_surf, &
32  alb, emis, cpAnOHM, kkAnOHM, chAnOHM, &! input
33  sfr, nsurf, EmissionsMethod, id, Gridiv, &
34  qn1_av_next, dqndt_next, &
35  a1, a2, a3, qs, deltaQi)! output
36 
37  IMPLICIT NONE
38  INTEGER, INTENT(in) :: tstep ! time step [s]
39  INTEGER, INTENT(in) :: dt_since_start ! time since simulation starts [s]
40 
41  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :)::MetForcingData_grid
42 
43  REAL(KIND(1d0)), INTENT(in):: qn1
44  REAL(KIND(1d0)), INTENT(in):: qf
45  REAL(KIND(1d0)), INTENT(in):: sfr(nsurf)
46  REAL(KIND(1d0)), INTENT(in):: moist_surf(nsurf)
47 
48  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::alb
49  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::emis
50  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::cpAnOHM
51  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::kkAnOHM
52  REAL(KIND(1d0)), INTENT(in), DIMENSION(:)::chAnOHM
53 
54  INTEGER, INTENT(in):: id
55  INTEGER, INTENT(in):: Gridiv
56  INTEGER, INTENT(in):: EmissionsMethod
57  INTEGER, INTENT(in):: nsurf
58  ! INTEGER,INTENT(in):: nsh !< number of timesteps in one hour [-]
59 
60  REAL(KIND(1d0)), INTENT(in)::qn1_av_prev
61  REAL(KIND(1d0)), INTENT(out)::qn1_av_next
62  REAL(KIND(1d0)), INTENT(in)::dqndt_prev !Rate of change of net radiation [W m-2 h-1] at t-1
63  REAL(KIND(1d0)), INTENT(out)::dqndt_next !Rate of change of net radiation [W m-2 h-1] at t-1
64 
65  ! REAL(KIND(1d0)),INTENT(inout)::qn1_store(nsh) !< stored qn1 [W m-2]
66  ! REAL(KIND(1d0)),INTENT(inout)::qn1_av_store(2*nsh+1) !< average net radiation over previous hour [W m-2]
67 
68  REAL(KIND(1d0)), INTENT(out):: a1
69  REAL(KIND(1d0)), INTENT(out):: a2
70  REAL(KIND(1d0)), INTENT(out):: a3
71  REAL(KIND(1d0)), INTENT(out):: qs
72  REAL(KIND(1d0)), INTENT(out):: deltaQi(nsurf)
73 
74  INTEGER :: is, xid
75  INTEGER, SAVE :: id_save ! store index of the valid day with enough data
76  REAL(KIND(1d0)), PARAMETER::NotUsed = -55.5
77  INTEGER, PARAMETER::notUsedI = -55
78  LOGICAL :: idQ ! whether id contains enough data
79 
80  ! REAL(KIND(1d0)) :: dqndt !< rate of change of net radiation [W m-2 h-1] at t-2
81  ! REAL(KIND(1d0)) :: surfrac !< surface fraction accounting for SnowFrac if appropriate
82  REAL(KIND(1d0)), DIMENSION(nsurf) :: xa1, xa2, xa3
83  ! REAL(KIND(1d0)) :: qn1_av ! average net radiation over previous hour [W m-2]
84  ! REAL(KIND(1d0)) :: nsh_nna ! number of timesteps per hour with non -999 values (used for spinup)
85 
86  ! initialize output variables
87  xa1 = 0.1
88  xa2 = 0.2
89  xa3 = 10
90  qs = -999
91  deltaqi = 0 ! NB: as snow part is not implemented within AnOHM yet, this is just a placeholder
92 
93  ! to test if the current met block contains enough data for AnOHM
94  ! TODO: more robust selection should be implemented
95  ! daylight hours >= 6
96  idq = count(metforcingdata_grid(:, 2) == id .AND. & ! day of year
97  metforcingdata_grid(:, 4) == 0 .AND. & ! minutes
98  metforcingdata_grid(:, 15) > 0) & ! Sd
99  >= 6
100 
101  ! PRINT*, idQ
102  IF (idq) THEN
103  ! given enough data, calculate coefficients of day `id`
104  xid = id
105  id_save = id ! store index of the valid day with enough data
106  ELSE
107  ! otherwise calculate coefficients of yesterday: `id-1`
108  xid = id_save
109  END IF
110 
111  DO is = 1, nsurf
112  ! call AnOHM to calculate the coefs.
113  CALL anohm_coef(is, xid, gridiv, metforcingdata_grid, moist_surf, emissionsmethod, qf, & !input
114  alb, emis, cpanohm, kkanohm, chanohm, &! input
115  xa1(is), xa2(is), xa3(is)) ! output
116  ! print*, 'AnOHM_coef are: ',xa1,xa2,xa3
117  ENDDO
118 
119  ! calculate the areally-weighted OHM coefficients
120  a1 = dot_product(xa1, sfr)
121  a2 = dot_product(xa2, sfr)
122  a3 = dot_product(xa3, sfr)
123 
124  ! Calculate radiation part ------------------------------------------------------------
125  qs = -999 !qs = Net storage heat flux [W m-2]
126  IF (qn1 > -999) THEN !qn1 = Net all-wave radiation [W m-2]
127 
128  ! Store instantaneous qn1 values for previous hour (qn1_store) and average (qn1_av)
129  ! CALL OHM_dqndt_cal(nsh,qn1,qn1_store,qn1_av_store,dqndt)
130  CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn1_av_prev, qn1, dqndt_prev, &
131  qn1_av_next, dqndt_next)
132 
133  ! Calculate net storage heat flux
134  CALL ohm_qs_cal(qn1, dqndt_prev, a1, a2, a3, qs)
135 
136  ELSE
137  CALL errorhint(21, 'SUEWS_AnOHM.f95: bad value for qn found during qs calculation.', qn1, notused, notusedi)
138  ENDIF
139 
140  END SUBROUTINE anohm
141  !========================================================================================
142 
143  !========================================================================================
149  SUBROUTINE anohm_coef( &
150  sfc_typ, xid, xgrid, &!input
151  MetForcingData_grid, moist, EmissionsMethod, qf, & !input
152  alb, emis, cpAnOHM, kkAnOHM, chAnOHM, &! input
153  xa1, xa2, xa3) ! output
155  IMPLICIT NONE
156 
157  ! input
158  INTEGER, INTENT(in):: sfc_typ
159  INTEGER, INTENT(in):: xid
160  INTEGER, INTENT(in):: xgrid
161  INTEGER, INTENT(in):: EmissionsMethod
162 
163  REAL(KIND(1d0)), INTENT(in):: qf
164  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: alb
165  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: emis
166  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: cpAnOHM
167  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: kkAnOHM
168  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: chAnOHM
169  REAL(KIND(1d0)), INTENT(in), DIMENSION(:) :: moist
170  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :) :: MetForcingData_grid
171 
172  ! output
173  REAL(KIND(1d0)), INTENT(out) :: xa1
174  REAL(KIND(1d0)), INTENT(out) :: xa2
175  REAL(KIND(1d0)), INTENT(out) :: xa3
176 
177  ! surface temperature related scales:
178  REAL(KIND(1d0)):: ATs
179  REAL(KIND(1d0)):: mTs
180  REAL(KIND(1d0)):: gamma
181 
182  ! forcing scales
183  REAL(KIND(1d0))::ASd, mSd, tSd
184  REAL(KIND(1d0))::ATa, mTa, tTa
185  REAL(KIND(1d0))::tau
186  REAL(KIND(1d0))::mWS, mWF, mAH
187 
188  ! forcings:
189  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::Sd ! incoming solar radiation [W m-2]
190  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::Ta ! air temperature [degC]
191  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::RH ! relative humidity [%]
192  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::pres ! air pressure [hPa]
193  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::WS ! wind speed [m s-1]
194  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::WF ! water flux density [m3 s-1 m-2]
195  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::AH ! anthropogenic heat [W m-2]
196  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE::tHr ! time [hr]
197 
198  ! sfc. properties:
199  REAL(KIND(1d0)) ::xalb ! albedo,
200  REAL(KIND(1d0)) ::xemis ! emissivity,
201  REAL(KIND(1d0)) ::xcp ! heat capacity,
202  REAL(KIND(1d0)) ::xk ! thermal conductivity,
203  REAL(KIND(1d0)) ::xch ! bulk transfer coef.
204  REAL(KIND(1d0)) ::xBo ! Bowen ratio
205  REAL(KIND(1d0)) ::xeta ! effective absorption coefficient
206  REAL(KIND(1d0)) ::xmu ! effective absorption fraction
207  REAL(KIND(1d0)) ::xmoist ! surface wetness
208 
209  ! locally saved variables:
210  ! if coefficients have been calculated, just reload them
211  ! otherwise, do the calculation
212  INTEGER, SAVE :: id_save, grid_save
213  REAL(KIND(1d0)), SAVE:: coeff_grid_day(7, 3) = -999.
214 
215  ! PRINT*, 'xid,id_save',xid,id_save
216  ! PRINT*, 'xgrid,grid_save',xgrid,grid_save
217  ! PRINT*, 'sfc_typ',sfc_typ
218  ! PRINT*, 'coeff_grid_day',coeff_grid_day(sfc_typ,:)
219  IF (xid == id_save .AND. xgrid == grid_save) THEN
220  ! if coefficients have been calculated, just reload them
221  ! print*, 'here no repetition'
222  xa1 = coeff_grid_day(sfc_typ, 1)
223  xa2 = coeff_grid_day(sfc_typ, 2)
224  xa3 = coeff_grid_day(sfc_typ, 3)
225  ELSE
226 
227  ! load forcing characteristics:
228  CALL anohm_fc( &
229  xid, metforcingdata_grid, emissionsmethod, qf, & ! input
230  asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah) ! output
231 
232  ! load forcing variables:
233  CALL anohm_fcload( &
234  xid, metforcingdata_grid, emissionsmethod, qf, & ! input
235  sd, ta, rh, pres, ws, wf, ah, thr) ! output
236 
237  ! load sfc. properties:
238  xalb = alb(sfc_typ)
239  xemis = emis(sfc_typ)
240  xcp = cpanohm(sfc_typ)
241  xk = kkanohm(sfc_typ)
242  xch = chanohm(sfc_typ)
243  xmoist = moist(sfc_typ)
244 
245  ! PRINT*, 'xBo before:',xBo
246  ! calculate Bowen ratio:
247  CALL anohm_bo_cal( &
248  sfc_typ, &
249  sd, ta, rh, pres, thr, & ! input: forcing
250  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
251  xalb, xemis, xcp, xk, xch, xmoist, & ! input: sfc properties
252  tsd, & ! input: peaking time of Sd in hour
253  xbo) ! output: Bowen ratio
254  ! PRINT*, 'xBo after:',xBo
255 
256  ! calculate AnOHM coefficients
257  SELECT CASE (sfc_typ)
258  CASE (1:6) ! land surfaces
259  CALL anohm_coef_land_cal( &
260  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
261  xalb, xemis, xcp, xk, xch, xbo, & ! input: sfc properties
262  xa1, xa2, xa3, ats, mts, gamma) ! output: surface temperature related scales by AnOHM
263 
264  CASE (7) ! water surface
265  ! NB:give fixed values for the moment
266  xeta = 0.3
267  xmu = 0.2
268  CALL anohm_coef_water_cal( &
269  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
270  xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, & ! input: sfc properties
271  xa1, xa2, xa3, ats, mts, gamma) ! output
272 
273  ! save variables for marking status as done
274  id_save = xid
275  grid_save = xgrid
276 
277  END SELECT
278  ! save variables for reusing values of the same day
279  coeff_grid_day(sfc_typ, :) = (/xa1, xa2, xa3/)
280  END IF
281 
282  END SUBROUTINE anohm_coef
283  !========================================================================================
284 
285  !========================================================================================
290  SUBROUTINE anohm_xts( &
291  sfc_typ, & !input: surface type
292  ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
293  xalb, xemis, xcp, xk, xch, xBo, & ! input: sfc properties
294  tSd, & !input: peaking time of Sd in hour
295  xTHr, &! input: time in hour
296  xTs)! output: surface temperature
298  IMPLICIT NONE
299  ! input:
300  INTEGER, INTENT(in):: sfc_typ
301 
302  ! input: forcing scales
303  REAL(KIND(1d0)), INTENT(in):: ASd
304  REAL(KIND(1d0)), INTENT(in):: mSd
305  REAL(KIND(1d0)), INTENT(in):: ATa
306  REAL(KIND(1d0)), INTENT(in):: mTa
307  REAL(KIND(1d0)), INTENT(in):: tau
308  REAL(KIND(1d0)), INTENT(in):: mWS
309  REAL(KIND(1d0)), INTENT(in):: mWF
310  REAL(KIND(1d0)), INTENT(in):: mAH
311 
312  ! input: sfc properties
313  REAL(KIND(1d0)), INTENT(in):: xalb
314  REAL(KIND(1d0)), INTENT(in):: xemis
315  REAL(KIND(1d0)), INTENT(in):: xcp
316  REAL(KIND(1d0)), INTENT(in):: xk
317  REAL(KIND(1d0)), INTENT(in):: xch
318  REAL(KIND(1d0)), INTENT(in):: xBo
319  REAL(KIND(1d0)):: xeta
320  REAL(KIND(1d0)):: xmu
321 
322  ! input: temporal-related
323  REAL(KIND(1d0)), INTENT(in):: tSd
324  REAL(KIND(1d0)), INTENT(in):: xTHr
325 
326  ! output:
327  REAL(KIND(1d0)), INTENT(out) :: xTs
328 
329  ! local
330  REAL(KIND(1d0)) :: &
331  xa1, xa2, xa3, &!coefficients
332  ATs, mTs, gamma !surface temperature related scales by AnOHM
333 
334  ! constant:
335  REAL(KIND(1d0)), PARAMETER :: PI = atan(1.0)*4 ! Pi
336  REAL(KIND(1d0)), PARAMETER :: OMEGA = 2*pi/(24*60*60) ! augular velocity of Earth
337 
338  SELECT CASE (sfc_typ)
339  CASE (1:6) ! land surfaces
340  CALL anohm_coef_land_cal( &
341  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
342  xalb, xemis, xcp, xk, xch, xbo, & ! input: sfc properties
343  xa1, xa2, xa3, ats, mts, gamma) ! output: surface temperature related scales by AnOHM
344 
345  CASE (7) ! water surface
346  ! ! NB:give fixed values for the moment
347  xeta = 0.3
348  xmu = 0.2
349  CALL anohm_coef_water_cal( &
350  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
351  xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, & ! input: sfc properties
352  xa1, xa2, xa3, ats, mts, gamma) ! output
353 
354  END SELECT
355 
356  ! for a local time xTHr (in hour):
357  xts = ats*sin(omega*(xthr - tsd + 6)*3600 - gamma) + mts
358 
359  END SUBROUTINE anohm_xts
360  !========================================================================================
361 
362  !========================================================================================
363  SUBROUTINE anohm_coef_land_cal( &
364  ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
365  xalb, xemis, xcp, xk, xch, xBo, & ! input: sfc properties
366  xa1, xa2, xa3, ATs, mTs, gamma) ! output
368  IMPLICIT NONE
369 
370  ! input: forcing scales
371  REAL(KIND(1d0)), INTENT(in):: ASd
372  REAL(KIND(1d0)), INTENT(in):: mSd
373  REAL(KIND(1d0)), INTENT(in):: ATa
374  REAL(KIND(1d0)), INTENT(in):: mTa
375  REAL(KIND(1d0)), INTENT(in):: tau
376  REAL(KIND(1d0)), INTENT(in):: mWS
377  REAL(KIND(1d0)), INTENT(in):: mWF
378  REAL(KIND(1d0)), INTENT(in):: mAH
379  ! input: sfc properties
380  REAL(KIND(1d0)), INTENT(in):: xalb
381  REAL(KIND(1d0)), INTENT(in):: xemis
382  REAL(KIND(1d0)), INTENT(in):: xcp
383  REAL(KIND(1d0)), INTENT(in):: xk
384  REAL(KIND(1d0)), INTENT(in):: xch
385  REAL(KIND(1d0)), INTENT(in):: xBo
386 
387  ! output
388  REAL(KIND(1d0)), INTENT(out) :: xa1
389  REAL(KIND(1d0)), INTENT(out) :: xa2
390  REAL(KIND(1d0)), INTENT(out) :: xa3
391  REAL(KIND(1d0)), INTENT(out) :: ATs
392  REAL(KIND(1d0)), INTENT(out) :: mTs
393  REAL(KIND(1d0)), INTENT(out) :: gamma
394 
395  ! constant
396  REAL(KIND(1d0)), PARAMETER :: SIGMA = 5.67e-8 ! Stefan-Boltzman
397  REAL(KIND(1d0)), PARAMETER :: PI = atan(1.0)*4 ! Pi
398  REAL(KIND(1d0)), PARAMETER :: OMEGA = 2*pi/(24*60*60) ! augular velocity of Earth
399 
400  ! local variables:
401  REAL(KIND(1d0)) :: beta ! inverse Bowen ratio
402  REAL(KIND(1d0)) :: f, fL, fT ! energy redistribution factors
403  REAL(KIND(1d0)) :: lambda ! thermal diffusivity
404  REAL(KIND(1d0)) :: delta, m, n ! water flux related variables
405  REAL(KIND(1d0)) :: ms, ns ! m, n related
406  REAL(KIND(1d0)) :: ceta, cphi ! phase related temporary variables
407  REAL(KIND(1d0)) :: eta, phi, xlag ! phase related temporary variables
408  REAL(KIND(1d0)) :: xx1, xx2, xx3, xchWS ! temporary use
409  ! LOGICAL :: flagGood = .TRUE. ! quality flag, T for good, F for bad
410 
411  ! give fixed values for test
412  ! properties
413  ! xalb = .2
414  ! xemis = .9
415  ! xcp = 1e6
416  ! xk = 1.2
417  ! xch = 4
418  ! xBo = 2
419  ! forcings
420  ! ASd = 400
421  ! mSd = 100
422  ! ATa = 23
423  ! mTa = 23+C2K
424  ! tau = PI/6
425  ! WS = 1
426  ! AH = 0
427  ! Wf = 0
428 
429  ! PRINT*, '********sfc_typ: ',sfc_typ,' start********'
430  ! initialize flagGood as .TRUE.
431  ! flagGood = .TRUE.
432  ! PRINT*, flagGood
433 
434  ! calculate sfc properties related parameters:
435  xchws = xch*mws
436 
437  ! xch = CRA/RA !update bulk trsf. coeff. with RA (aerodyn. res.)
438  IF (abs(xbo) < 0.001) THEN
439  beta = 1/0.001
440  ELSE
441  beta = 1/xbo
442  END IF
443 
444  ! PRINT*, 'beta:', beta
445  f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
446  ! print*, 'xch',xch,'mTa',mTa,'dLu',4*SIGMA*xemis*mTa**3
447  fl = 4*sigma*xemis*mta**3
448  ! print*, 'fL',fL
449  ft = (1 + beta)*xchws
450  ! print*, 'fT',fT
451  lambda = xk/xcp
452  ! print*, 'lambda',lambda
453  delta = sqrt(.5*(mwf**2 + sqrt(mwf**4 + 16*lambda**2*omega**2)))
454  ! print*, 'delta',delta
455  m = (2*lambda)/(delta + mwf)
456  n = delta/omega
457  ! print*, 'm',m,'n',n
458 
459  ! calculate surface temperature related parameters:
460  mts = (msd*(1 - xalb)/f) + mta
461  ms = 1 + xk/(f*m)
462  ns = xk/(f*n)
463  ! print*, 'ms,ns:', ms,ns
464  xx1 = f*sin(tau)*ata
465  ! print*, 'xx1',xx1
466  xx2 = (1 - xalb)*asd + f*cos(tau)*ata
467  ! print*, 'xx2',xx2
468  gamma = atan(ns/ms) + atan(xx1/xx2)
469  ! print*, 'gamma:', gamma
470  ats = -(sin(tau)*ata)/(ns*cos(gamma) - ms*sin(gamma))
471  ! print*, 'ATs:', ATs
472 
473  ! calculate net radiation related parameters:
474  xx1 = (ns*cos(gamma) + sin(gamma) - ms*sin(gamma))*sin(tau)*ata*fl
475  xx2 = (xalb - 1)*(ns*cos(gamma) - ms*sin(gamma))*asd
476  xx3 = (-ms*cos(tau)*sin(tau) + cos(gamma)*(ns*cos(tau) + sin(tau)))*ata*fl
477  xx2 = xx2 - xx3
478  phi = atan(xx1/xx2)
479  xx3 = (ns*cos(gamma) - ms*sin(gamma))
480  xx1 = (1 + sin(gamma)/xx3)*sin(tau)*ata*fl
481  xx2 = (xalb - 1)*asd - (cos(tau) + cos(gamma)*sin(tau)/xx3)*ata*fl
482  cphi = sqrt(xx1**2 + xx2**2)
483 
484  ! calculate heat storage related parameters:
485  xx1 = m*cos(gamma) - n*sin(gamma)
486  xx2 = m*sin(gamma) + n*cos(gamma)
487  eta = atan(xx1/xx2)
488  ! if ( eta<0 ) print*, 'lambda,delta,m,n,gamma:', lambda,delta,m,n,gamma
489  xx1 = xk**2*(m**2 + n**2)*ats**2
490  xx2 = m**2*n**2
491  ceta = sqrt(xx1/xx2)
492 
493  ! key phase lag:
494  xlag = eta - phi
495 
496  ! calculate the OHM coeffs.:
497  ! a1:
498  xa1 = (ceta*cos(xlag))/cphi
499 
500  ! a2:
501  xa2 = (ceta*sin(xlag))/(omega*cphi)
502  xa2 = xa2/3600 ! convert the unit from s-1 to h-1
503 
504  ! a3:
505  xa3 = -xa1*(ft/f)*(msd*(1 - xalb)) - mah
506 
507  ! print*, 'ceta,cphi', ceta,cphi
508  ! print*, 'tau,eta,phi,xlag in deg:',tau/pi*180,eta/pi*180,phi/pi*180,xlag/pi*180
509 
510  ! PRINT*, '********sfc_typ: ',sfc_typ,' end********'
511 
512  END SUBROUTINE anohm_coef_land_cal
513  !========================================================================================
514 
515  !========================================================================================
520  SUBROUTINE anohm_coef_water_cal( &
521  ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
522  xalb, xemis, xcp, xk, xch, xBo, xeta, xmu, & ! input: sfc properties
523  xa1, xa2, xa3, ATs, mTs, gamma) ! output
525  IMPLICIT NONE
526 
527  ! input: forcing scales
528  REAL(KIND(1d0)), INTENT(in):: ASd
529  REAL(KIND(1d0)), INTENT(in):: mSd
530  REAL(KIND(1d0)), INTENT(in):: ATa
531  REAL(KIND(1d0)), INTENT(in):: mTa
532  REAL(KIND(1d0)), INTENT(in):: tau
533  REAL(KIND(1d0)), INTENT(in):: mWS
534  REAL(KIND(1d0)), INTENT(in):: mWF
535  REAL(KIND(1d0)), INTENT(in):: mAH
536 
537  ! input: sfc properties
538  REAL(KIND(1d0)), INTENT(in):: xalb
539  REAL(KIND(1d0)), INTENT(in):: xemis
540  REAL(KIND(1d0)), INTENT(in):: xcp
541  REAL(KIND(1d0)), INTENT(in):: xk
542  REAL(KIND(1d0)), INTENT(in):: xch
543  REAL(KIND(1d0)), INTENT(in):: xBo
544  REAL(KIND(1d0)), INTENT(in):: xeta
545  REAL(KIND(1d0)), INTENT(in):: xmu
546 
547  ! output
548  REAL(KIND(1d0)), INTENT(out) :: xa1
549  REAL(KIND(1d0)), INTENT(out) :: xa2
550  REAL(KIND(1d0)), INTENT(out) :: xa3
551  REAL(KIND(1d0)), INTENT(out) :: ATs
552  REAL(KIND(1d0)), INTENT(out) :: mTs
553  REAL(KIND(1d0)), INTENT(out) :: gamma
554 
555  ! constant
556  REAL(KIND(1d0)), PARAMETER :: SIGMA = 5.67e-8 ! Stefan-Boltzman
557  REAL(KIND(1d0)), PARAMETER :: PI = atan(1.0)*4 ! Pi
558  REAL(KIND(1d0)), PARAMETER :: OMEGA = 2*pi/(24*60*60) ! augular velocity of Earth
559 
560  ! local variables:
561  REAL(KIND(1d0)) :: beta ! inverse Bowen ratio
562  REAL(KIND(1d0)) :: f, fL, fT ! energy redistribution factors
563  REAL(KIND(1d0)) :: lambda, calb ! temporary use
564  REAL(KIND(1d0)) :: delta ! sfc. temperature related variables
565  ! REAL(KIND(1d0)) :: delta,m,n ! sfc. temperature related variables
566  REAL(KIND(1d0)) :: xm, xn ! m, n related
567  ! REAL(KIND(1d0)) :: gamma ! phase lag scale
568  REAL(KIND(1d0)) :: phi ! phase lag scale
569  ! REAL(KIND(1d0)) :: ATs,mTs ! surface temperature amplitude
570  REAL(KIND(1d0)) :: czeta, ctheta ! phase related temporary variables
571  REAL(KIND(1d0)) :: zeta, theta, xlag ! phase related temporary variables
572  REAL(KIND(1d0)) :: xx1, xx2, xx3 ! temporary use
573  REAL(KIND(1d0)) :: kappa ! temporary use
574  REAL(KIND(1d0)) :: dtau, dpsi, dphi ! temporary use
575  REAL(KIND(1d0)) :: cdtau, cdpsi, cdphi ! temporary use
576  REAL(KIND(1d0)) :: xxT, xxkappa, xxdltphi, xchWS ! temporary use
577  ! LOGICAL :: flagGood = .TRUE. ! quality flag, T for good, F for bad
578 
579  ! ====not used====
580  REAL(KIND(1d0)) :: dummy
581  dummy = mah + mwf
582  ! ====not used====
583 
584  ! calculate sfc properties related parameters:
585  xm = xk*xmu**2
586  xn = xcp*omega
587  phi = atan(xn/xm)
588  kappa = sqrt(xcp*omega/(2*xk))
589 
590  ! mWS=SUM(WS, dim=1, mask=(WS>0))/SIZE(WS, dim=1)
591  xchws = xch*mws
592  beta = 1/xbo
593  f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
594  fl = 4*sigma*xemis*mta**3
595  ft = (1 + beta)*xchws
596 
597  calb = 1 - xalb
598 
599  lambda = sqrt(xm**2 + xn**2)
600 
601  dtau = atan(xk*kappa/(f + xk*kappa))
602  dpsi = atan((xk - xmu)/kappa)
603  dphi = atan(kappa*(f + xk*xmu)/(f*(kappa - xmu) + xk*kappa*(2*kappa - xmu)))
604  cdtau = sqrt((xk*kappa)**2 + (f + xk*kappa)**2)
605  cdpsi = sqrt((xk - xmu)**2 + kappa**2)
606  cdphi = cdtau*cdpsi
607 
608  ! calculate surface temperature related parameters:
609  ! daily mean:
610  mts = (msd*(1 - xalb + xeta)/f) + mta
611  ! amplitude:
612  xx1 = (xk*xeta*xmu*calb*asd*cdpsi)**2
613  xx2 = 2*lambda*sqrt(xx1)*(calb*asd*sin(phi - dpsi) + f*ata*sin(tau + phi - dpsi))
614  xx3 = lambda**2*((calb*asd + cos(tau)*f*ata)**2 + (sin(tau)*f*ata)**2)
615  ats = 1/(cdtau*lambda)*sqrt(xx1 + xx2 + xx3)
616  ! phase lag:
617  xx1 = (xk*kappa*calb*asd + cdtau*f*ata*sin(tau + dtau))*lambda &
618  + (xk*xeta*xmu*calb*asd*cdphi*sin(phi + dphi))
619  xx2 = ((f + xk*kappa)*calb*asd - cdtau*f*ata*cos(tau + dtau))*lambda &
620  - xk*xeta*xmu*calb*asd*cdphi*cos(phi + dphi)
621  delta = atan(xx1/xx2)
622  gamma = delta
623 
624  ! calculate net radiation related parameters:
625  ! phase lag:
626  xx1 = fl*(ats*sin(delta) - ata*sin(tau))
627  xx2 = calb*asd - fl*(ats*cos(delta) - ata*cos(tau))
628  theta = atan(xx1/xx2)
629  ! amplitude:
630  ctheta = sqrt(xx1**2 + xx2**2)
631 
632  ! calculate heat storage related parameters:
633  ! scales:
634  xxt = sqrt(2.)*kappa*lambda*ats
635  xxkappa = cdpsi*xeta*xmu*asd
636  xxdltphi = cos(delta)*sin(dpsi)*cos(phi) - sin(delta)*cos(dpsi)*sin(phi)
637  ! phase lag:
638  xx1 = xxt*sin(pi/4 - delta) + xxkappa*sin(phi + dpsi)
639  xx2 = xxt*sin(pi/4 + delta) - xxkappa*sin(phi - dpsi)
640  zeta = atan(xx1/xx2)
641  ! amplitude:
642  xx1 = 2*sqrt(2.)*xxkappa*xxt*xxdltphi
643  xx2 = (1 - cos(2*dpsi)*cos(2*phi))*xxkappa**2
644  xx3 = xxt**2
645  czeta = xk/lambda*sqrt(xx1 + xx2 + xx3)
646 
647  ! calculate the OHM coeffs.:
648  xlag = zeta - theta
649  ! a1:
650  xa1 = (czeta*cos(xlag))/ctheta
651 
652  ! write(*,*) 'ceta,xlag,cphi:', ceta,xlag,cphi
653  ! a2:
654  xa2 = (czeta*sin(xlag))/(omega*ctheta)
655  xa2 = xa2/3600 ! convert the unit from s-1 to h-1
656 
657  ! a3:
658  xa3 = msd*(xalb - 1)*(xeta + (ft - fl*xeta)/f*xa1)
659 
660  ! print*, 'sfc_typ_water:', sfc_typ
661  ! print*, 'a1,a2,a3:', xa1,xa2,xa3
662 
663  END SUBROUTINE anohm_coef_water_cal
664  !========================================================================================
665 
666  !========================================================================================
667  SUBROUTINE anohm_fc( &
668  xid, MetForcingData_grid, EmissionsMethod, qf, & ! input
669  ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH) ! output
671  IMPLICIT NONE
672 
673  ! input
674  INTEGER, INTENT(in):: xid
675  INTEGER, INTENT(in):: EmissionsMethod
676  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :) ::MetForcingData_grid
677  REAL(KIND(1d0)), INTENT(in):: qf
678 
679  ! output
680  REAL(KIND(1d0)), INTENT(out):: ASd
681  REAL(KIND(1d0)), INTENT(out):: mSd
682  REAL(KIND(1d0)), INTENT(out):: tSd
683  REAL(KIND(1d0)), INTENT(out):: ATa
684  REAL(KIND(1d0)), INTENT(out):: mTa
685  REAL(KIND(1d0)), INTENT(out):: tTa
686  REAL(KIND(1d0)), INTENT(out):: tau
687  REAL(KIND(1d0)), INTENT(out):: mWS
688  REAL(KIND(1d0)), INTENT(out):: mWF
689  REAL(KIND(1d0)), INTENT(out):: mAH
690 
691  ! forcings:
692  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: Sd
693  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: Ta
694  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: RH
695  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: pres
696  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: WS ! wind speed [m s-1]
697  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: WF ! water flux density [m3 s-1 m-2]
698  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: AH ! anthropogenic heat [W m-2]
699  REAL(KIND(1d0)), DIMENSION(:), ALLOCATABLE:: tHr ! local time [hr]
700 
701  ! load forcing variables:
702  CALL anohm_fcload(xid, metforcingdata_grid, emissionsmethod, qf, sd, ta, rh, pres, ws, wf, ah, thr)
703  ! calculate forcing scales for AnOHM:
704  CALL anohm_fccal(sd, ta, ws, wf, ah, thr, asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
705 
706  ! CALL r8vec_print(SIZE(sd, dim=1),sd,'Sd')
707  ! PRINT*, ASd,mSd,tSd
708 
709  ! CALL r8vec_print(SIZE(ta, dim=1),ta,'Ta')
710  ! PRINT*, ATa,mTa,tTa
711 
712  END SUBROUTINE anohm_fc
713  !========================================================================================
714 
715  !========================================================================================
717  SUBROUTINE anohm_fcload( &
718  xid, MetForcingData_grid, EmissionsMethod, qf, & ! input
719  Sd, Ta, RH, pres, WS, WF, AH, tHr) ! output
721  IMPLICIT NONE
722 
723  ! input
724  INTEGER, INTENT(in):: xid
725  INTEGER, INTENT(in):: EmissionsMethod
726  REAL(KIND(1d0)), INTENT(in), DIMENSION(:, :) ::MetForcingData_grid
727  REAL(KIND(1d0)), INTENT(in):: qf
728 
729  ! output
730  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: Sd
731  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: Ta
732  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: RH
733  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: pres
734  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: WS
735  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: WF
736  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: AH
737  REAL(KIND(1d0)), DIMENSION(:), INTENT(out), ALLOCATABLE:: tHr
738 
739  ! local variables:
740  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE :: subMet ! subset array of daytime series
741 
742  INTEGER :: err
743  INTEGER :: lenMetData, nVar
744 
745  LOGICAL, ALLOCATABLE :: metMask(:)
746 
747  ! construct mask
748  IF (ALLOCATED(metmask)) DEALLOCATE (metmask, stat=err)
749  ALLOCATE (metmask(SIZE(metforcingdata_grid, dim=1)))
750  metmask = (metforcingdata_grid(:, 2) == xid & ! day=xid
751  .AND. metforcingdata_grid(:, 4) == 0)! tmin=0
752 
753  ! determine the length of subset
754  lenmetdata = count(metmask)
755 
756  ! construct array for time and met variables
757  nvar = 8! number of variables to retrieve
758  ! print*, 'good 1'
759  ! allocate subMet:
760  IF (ALLOCATED(submet)) DEALLOCATE (submet, stat=err)
761  ALLOCATE (submet(lenmetdata, nvar))
762  submet = reshape(pack(metforcingdata_grid(:, (/3, & !time: hour
763  15, 12, 11, 13, 10, 12, 9/)), &! met: Sd, Ta, RH, pres, WS, WF, AH
764  ! NB: WF not used: a placeholder
765  spread(metmask, dim=2, ncopies=nvar)), & ! replicate mask vector to 2D array
766  (/lenmetdata, nvar/)) ! convert to target shape
767 
768  ! re-allocate arrays as their sizes may change during passing
769  IF (ALLOCATED(thr)) DEALLOCATE (thr, stat=err)
770  ALLOCATE (thr(lenmetdata))
771  IF (ALLOCATED(sd)) DEALLOCATE (sd, stat=err)
772  ALLOCATE (sd(lenmetdata))
773  IF (ALLOCATED(ta)) DEALLOCATE (ta, stat=err)
774  ALLOCATE (ta(lenmetdata))
775  IF (ALLOCATED(rh)) DEALLOCATE (rh, stat=err)
776  ALLOCATE (rh(lenmetdata))
777  IF (ALLOCATED(pres)) DEALLOCATE (pres, stat=err)
778  ALLOCATE (pres(lenmetdata))
779  IF (ALLOCATED(ws)) DEALLOCATE (ws, stat=err)
780  ALLOCATE (ws(lenmetdata))
781  IF (ALLOCATED(wf)) DEALLOCATE (wf, stat=err)
782  ALLOCATE (wf(lenmetdata))
783  IF (ALLOCATED(ah)) DEALLOCATE (ah, stat=err)
784  ALLOCATE (ah(lenmetdata))
785 
786  ! load the sublist into forcing variables
787  thr = submet(:, 1)! time in hour
788  sd = submet(:, 2)
789  ta = submet(:, 3)
790  rh = submet(:, 4)
791  pres = submet(:, 5)
792  ws = submet(:, 6)
793  wf = 0 ! set as 0 for the moment
794  IF (emissionsmethod == 0) THEN
795  ah = submet(:, 8) ! read in from MetForcingData_grid,
796  ELSE
797  ! AH = 0 ! temporarily change to zero;
798  ah = qf ! use modelled value
799  ! AH = mAH_grids(xid-1,xgrid)
800  END IF
801 
802  END SUBROUTINE anohm_fcload
803  !========================================================================================
804 
805  !========================================================================================
808  SUBROUTINE anohm_fccal( &
809  Sd, Ta, WS, WF, AH, tHr, & ! input
810  ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH) ! output
811  IMPLICIT NONE
812 
813  ! input
814  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: Sd
815  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: Ta
816  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: WS
817  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: WF
818  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: AH
819  REAL(KIND(1d0)), DIMENSION(:), INTENT(in):: tHr
820 
821  ! output
822  REAL(KIND(1d0)), INTENT(out):: ASd
823  REAL(KIND(1d0)), INTENT(out):: mSd
824  REAL(KIND(1d0)), INTENT(out):: tSd
825  REAL(KIND(1d0)), INTENT(out):: ATa
826  REAL(KIND(1d0)), INTENT(out):: mTa
827  REAL(KIND(1d0)), INTENT(out):: tTa
828  REAL(KIND(1d0)), INTENT(out):: tau
829  REAL(KIND(1d0)), INTENT(out):: mWS
830  REAL(KIND(1d0)), INTENT(out):: mWF
831  REAL(KIND(1d0)), INTENT(out):: mAH
832 
833  ! constant
834  REAL(KIND(1d0)), PARAMETER :: PI = atan(1.0)*4 ! Pi
835  REAL(KIND(1d0)), PARAMETER :: C2K = 273.15 ! degC to K
836 
837  ! local variables:
838  REAL(KIND(1d0)), ALLOCATABLE :: tHrDay(:) ! daytime tHr when Sd>0
839  REAL(KIND(1d0)), ALLOCATABLE :: selX(:) ! daytime sublist of met variable when Sd>0
840 
841  ! REAL(KIND(1d0)) :: xx ! temporary use
842  INTEGER :: err, lenDay
843  LOGICAL, DIMENSION(:), ALLOCATABLE::SdMask
844 
845  ALLOCATE (sdmask(SIZE(sd, dim=1)), stat=err)
846  IF (err /= 0) print *, "SdMask: Allocation request denied"
847  sdmask = sd > 5
848  lenday = count(sdmask)
849 
850  ! CALL r8vec_print(24,Sd,'Sd')
851  ! CALL r8vec_print(24,tHr,'tHr')
852  ALLOCATE (thrday(lenday), stat=err)
853  IF (err /= 0) print *, "tHrDay: Allocation request denied"
854  thrday = pack(thr, mask=sdmask)
855 
856  ALLOCATE (selx(lenday), stat=err)
857  IF (err /= 0) print *, "selX: Allocation request denied"
858 
859  ! calculate sinusoidal scales of Sd:
860  ! PRINT*, 'Calc. Sd...'
861  selx = pack(sd, mask=sdmask)
862  asd = (maxval(selx) - minval(selx))/2
863  msd = sum(selx)/lenday
864  tsd = 12
865  CALL anohm_shapefit(thrday, selx, asd, msd, tsd)
866  ! CALL r8vec_print(lenDay,selX,'Sd Day:')
867  ! modify ill-shaped days to go through
868  ! IF ( ASd < 0 .OR. tSd > 15) THEN
869  ! ! ASd = abs(ASd)
870  ! ! tSd = 12 ! assume Sd peaks at 12:00LST
871  ! CALL r8vec_print(lenDay,tHrDay,'tHrDay:')
872  ! CALL r8vec_print(lenDay,selX,'Sd Day:')
873  ! PRINT*, 'ASd:', ASd
874  ! PRINT*, 'mSd:', mSd
875  ! PRINT*, 'tSd:', tSd
876  ! END IF
877  ! PRINT*, 'Sd Day:', selX
878 
879  ! add corrections on mSd to fix the overestimated a3 during cold periods
880  IF (msd + asd < maxval(selx)) THEN
881  ! PRINT*, 'Sd max:', MAXVAL(selX)
882  ! PRINT*, 'ASd:', ASd
883  ! PRINT*, 'mSd before:', mSd
884  msd = maxval(selx) - asd
885  ! PRINT*, 'mSd after:', mSd
886  ! PRINT*,''
887  END IF
888 
889  ! calculate sinusoidal scales of Ta:
890  ! PRINT*, 'Calc. Ta...'
891  selx = pack(ta, mask=sdmask)
892  ata = (maxval(selx) - minval(selx))/2
893  mta = sum(selx)/lenday
894  tta = 12
895  CALL anohm_shapefit(thrday, selx, ata, mta, tta)
896  ! CALL r8vec_print(lenDay,selX,'Ta Day:')
897  IF (mta < 60) mta = mta + c2k ! correct the Celsius to Kelvin
898  ! modify ill-shaped days to go through
899  IF (ata < 0) THEN
900  ! ATa = abs(ATa)
901  ! tTa = 14 ! assume Ta peaks at 14:00LST
902  CALL r8vec_print(lenday, selx, 'Ta Day:')
903  print *, 'ATa:', ata
904  print *, 'mTa:', mta
905  print *, 'tTa:', tta
906  END IF
907  ! PRINT*, 'Ta:', Ta(10:16)
908 
909  ! calculate the phase lag between Sd and Ta:
910  tau = (tta - tsd)/24*2*pi
911  ! PRINT*, 'tau:', tau
912 
913  ! calculate the mean values:
914  selx = pack(ws, mask=sdmask)
915  mws = sum(selx)/lenday ! mean value of WS
916 
917  selx = pack(wf, mask=sdmask)
918  mwf = sum(selx)/lenday ! mean value of WF
919 
920  selx = pack(ah, mask=sdmask)
921  mah = sum(selx)/lenday ! mean value of AH
922  ! PRINT*, 'mWS:', mWS
923  ! print*, 'mWF:', mWF
924  ! print*, 'mAH:', mAH
925  IF (ALLOCATED(sdmask)) DEALLOCATE (sdmask, stat=err)
926  IF (err /= 0) print *, "SdMask: Deallocation request denied"
927 
928  IF (ALLOCATED(thrday)) DEALLOCATE (thrday, stat=err)
929  IF (err /= 0) print *, "tHrDay: Deallocation request denied"
930 
931  IF (ALLOCATED(selx)) DEALLOCATE (selx, stat=err)
932  IF (err /= 0) print *, "selX: Deallocation request denied"
933 
934  END SUBROUTINE anohm_fccal
935  !========================================================================================
936 
937  !========================================================================================
940  SUBROUTINE anohm_shapefit( &
941  tHr, obs, & !input
942  amp, mean, tpeak)!output
943  IMPLICIT NONE
944 
945  ! input
946  REAL(KIND(1d0)), INTENT(in) :: tHr(:)
947  REAL(KIND(1d0)), INTENT(in) :: obs(:)
948 
949  ! output
950  REAL(KIND(1d0)), INTENT(out) :: amp
951  REAL(KIND(1d0)), INTENT(out) :: mean
952  REAL(KIND(1d0)), INTENT(out) :: tpeak
953 
954  INTEGER :: m, n, info, err ! temporary use
955 
956  ! EXTERNAL fSin
957  REAL(KIND(1d0)), ALLOCATABLE:: fvec(:), x(:)
958 
959  REAL(KIND(1d0)):: tol = 0.00001d+00 ! tolerance
960 
961  n = 3 ! number of parameters to determine
962  m = SIZE(thr, dim=1) ! number of observation pairs
963  ! PRINT*, 'm',m
964 
965  ALLOCATE (fvec(m), stat=err)
966  IF (err /= 0) print *, "fvec: Allocation request denied"
967 
968  ALLOCATE (x(n), stat=err)
969  IF (err /= 0) print *, "x: Allocation request denied"
970 
971  ! initial guess for fitting
972  x = (/mean, amp, tpeak/)
973 
974  ! use minpack subroutine 'lmstr1' to fit the sinusoidal form
975  CALL lmdif1(fsin, m, n, x, thr, obs, fvec, tol, info)
976 
977  ! x = (/mean,amp,delta/), see subroutine 'fSin' for the meaning of this vector
978  mean = x(1)
979  amp = x(2)
980  tpeak = x(3) + 6 ! when t = delta + 6 (hr of LST), t is tpeak.
981 
982  ! adjust fitted parameters to physical range:
983  ! amp>0
984  IF (amp < 0) THEN
985  ! PRINT*, ''
986  ! PRINT*, 'before:'
987  ! PRINT*, amp,tpeak,mean
988  ! CALL r8vec_print(SIZE(tHR, dim=1),tHR,'tHR')
989  ! CALL r8vec_print(SIZE(obs, dim=1),obs,'obs')
990  amp = abs(amp)
991  tpeak = x(3) - 12 + 6 + 24
992  tpeak = mod(tpeak, 24.)
993  ! PRINT*, 'after:'
994  ! PRINT*, amp,tpeak,mean
995  END IF
996  tpeak = mod(tpeak, 24.)
997 
998  END SUBROUTINE anohm_shapefit
999  !========================================================================================
1000 
1001  !========================================================================================
1005  SUBROUTINE fsin(m, n, x, xdat, ydat, fvec, iflag)
1007  IMPLICIT NONE
1008  INTEGER(kind=4) m
1009  INTEGER(kind=4) n
1010 
1011  ! REAL ( kind = 8 ) fjrow(n)
1012  REAL(kind=8) fvec(m), & ! residual vector
1013  xdat(m), ydat(m) ! (x,y) observations for fitting
1014 
1015  INTEGER(kind=4) iflag, i
1016  REAL(kind=8) x(n)
1017  REAL(KIND(1d0)), PARAMETER :: PI = atan(1.0)*4 ! Pi
1018 
1019  IF (iflag == 0) THEN
1020 
1021  WRITE (*, '(a)') ''
1022  DO i = 1, n
1023  WRITE (*, '(g14.6)') x(i)
1024  END DO
1025 
1026  ELSE IF (iflag == 1) THEN
1027  fvec(1:m) = x(1) + x(2)*sin(2*pi/24*(xdat(1:m) - x(3))) - ydat(1:m)
1028 
1029  END IF
1030  RETURN
1031 
1032  END SUBROUTINE fsin
1033  !========================================================================================
1034 
1035  !========================================================================================
1037  SUBROUTINE anohm_bo_cal( &
1038  sfc_typ, & ! surface type
1039  Sd, Ta, RH, pres, tHr, & ! input: forcing
1040  ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
1041  xalb, xemis, xcp, xk, xch, xSM, & ! input: sfc properties
1042  tSd, & ! input: peaking time of Sd in hour
1043  xBo) ! output: Bowen ratio
1045  IMPLICIT NONE
1046 
1047  ! input:
1048  INTEGER, INTENT(in) :: sfc_typ ! unknown Bowen ratio
1049 
1050  ! input: daytime series
1051  REAL(kind=8), INTENT(in), DIMENSION(:) ::Sd
1052  REAL(kind=8), INTENT(in), DIMENSION(:) ::Ta
1053  REAL(kind=8), INTENT(in), DIMENSION(:) ::RH
1054  REAL(kind=8), INTENT(in), DIMENSION(:) ::pres
1055  REAL(kind=8), INTENT(in), DIMENSION(:) ::tHr
1056 
1057  ! input: forcing scales
1058  REAL(KIND(1d0)), INTENT(in):: ASd
1059  REAL(KIND(1d0)), INTENT(in):: mSd
1060  REAL(KIND(1d0)), INTENT(in):: tSd
1061  REAL(KIND(1d0)), INTENT(in):: ATa
1062  REAL(KIND(1d0)), INTENT(in):: mTa
1063  REAL(KIND(1d0)), INTENT(in):: tau
1064  REAL(KIND(1d0)), INTENT(in):: mWS
1065  REAL(KIND(1d0)), INTENT(in):: mWF
1066  REAL(KIND(1d0)), INTENT(in):: mAH
1067 
1068  ! input: surface properties
1069  REAL(KIND(1d0)), INTENT(in) :: xalb
1070  REAL(KIND(1d0)), INTENT(in) :: xemis
1071  REAL(KIND(1d0)), INTENT(in) :: xcp
1072  REAL(KIND(1d0)), INTENT(in) :: xk
1073  REAL(KIND(1d0)), INTENT(in) :: xch
1074  REAL(KIND(1d0)), INTENT(in) :: xSM
1075 
1076  ! output:
1077  REAL(kind=8), INTENT(out) :: xBo ! unknown Bowen ratio [-]
1078 
1079  ! EXTERNAL fcnBo
1080 
1081  REAL(kind=8), ALLOCATABLE :: x(:), fvec(:), prms(:)
1082 
1083  INTEGER :: lenDay, n, m, info, err, nVar, nPrm
1084  LOGICAL, DIMENSION(:), ALLOCATABLE :: maskDay
1085  REAL(kind=8) :: tol = 1e-20
1086 
1087  ! LOGICAL, ALLOCATABLE,DIMENSION(:) :: metMask
1088 
1089  ! daytime mask
1090  ALLOCATE (maskday(SIZE(sd)), stat=err)
1091  maskday = sd > 5
1092  ! length of daytime series
1093  lenday = SIZE(pack(sd, mask=maskday), dim=1)
1094  ! ! daytime mask
1095  ! IF (ALLOCATED(metMask)) DEALLOCATE(metMask, stat=err)
1096  ! ALLOCATE(metMask(lenDay))
1097  ! metMask=(Sd>0)
1098 
1099  ! assign x vector:
1100  n = 1
1101  ALLOCATE (x(n), stat=err)
1102  IF (err /= 0) print *, "x: Allocation request denied"
1103  ALLOCATE (fvec(n), stat=err)
1104  IF (err /= 0) print *, "fvec: Allocation request denied"
1105 
1106  ! pass initial Bowen ratio:
1107  xbo = 10
1108  x(1) = xbo
1109 
1110  ! NB: set these numbers properly if any changes made to this subroutine:
1111  ! number of parameters to pass:
1112  nprm = 16
1113  ! number of met variables to pass:
1114  nvar = 5
1115  ! length of x vector that holds unknown and parameters
1116  m = nprm + nvar*lenday
1117 
1118  ALLOCATE (prms(m), stat=err)
1119  IF (err /= 0) print *, "prms: Allocation request denied"
1120 
1121  ! pass forcing scales:
1122  prms(1) = asd
1123  prms(2) = msd
1124  prms(3) = ata
1125  prms(4) = mta
1126  prms(5) = tau
1127  prms(6) = mws
1128  prms(7) = mwf
1129  prms(8) = mah
1130 
1131  ! pass sfc. property scales:
1132  prms(9) = xalb
1133  prms(10) = xemis
1134  prms(11) = xcp
1135  prms(12) = xk
1136  prms(13) = xch
1137  prms(14) = xsm
1138 
1139  ! pass tSd:
1140  prms(15) = tsd
1141 
1142  ! pass tSd:
1143  prms(16) = sfc_typ*1.0
1144 
1145  ! extract daytime series
1146  prms(nprm + 1:m) = pack((/sd, ta, rh, pres, thr/), &
1147  mask=pack(spread(maskday, dim=2, ncopies=nvar), .true.))
1148 
1149  ! PRINT*, 'xBo before solve:',x(1)
1150  ! PRINT*, 'fvec before solve:',fvec(1)
1151  ! PRINT*, 'xSM:',xSM
1152  ! solve nonlinear equation fcnBo(x)=0
1153  CALL hybrd1(fcnbo, n, x, fvec, tol, info, m, prms)
1154  xbo = x(1)
1155  ! PRINT*, 'xBo after solve: ',x(1)
1156  ! PRINT*, 'fvec after solve:',fvec(1)
1157 
1158  IF (ALLOCATED(x)) DEALLOCATE (x, stat=err)
1159  IF (err /= 0) print *, "x: Deallocation request denied"
1160  IF (ALLOCATED(fvec)) DEALLOCATE (fvec, stat=err)
1161  IF (err /= 0) print *, "fvec: Deallocation request denied"
1162  IF (ALLOCATED(prms)) DEALLOCATE (prms, stat=err)
1163  IF (err /= 0) print *, "prms: Deallocation request denied"
1164 
1165  END SUBROUTINE anohm_bo_cal
1166  !========================================================================================
1167 
1168  !========================================================================================
1170  SUBROUTINE fcnbo(n, x, fvec, iflag, m, prms)
1172  IMPLICIT NONE
1173  ! Input, external FCN, the name of the user-supplied subroutine which
1174  ! calculates the functions. The routine should have the form:
1175  ! subroutine fcn ( n, x, fvec, iflag )
1176  INTEGER(kind=4):: n
1177  INTEGER(kind=4):: m
1178  INTEGER(kind=4):: iflag
1179  REAL(kind=8):: fvec(n)
1180  REAL(kind=8):: x(n) ! x(1) as unknown Bo
1181  REAL(kind=8):: prms(m) ! prms(i) used for passing parameters
1182 
1183  ! the unknow: Bowen ratio
1184  REAL(kind=8):: xBo
1185 
1186  ! forcing scales
1187  REAL(kind=8):: ASd
1188  REAL(kind=8):: mSd
1189  REAL(kind=8):: ATa
1190  REAL(kind=8):: mTa
1191  REAL(kind=8):: tau
1192  REAL(kind=8):: mWS
1193  REAL(kind=8):: mWF
1194  REAL(kind=8):: mAH
1195 
1196  ! sfc. property scales
1197  REAL(kind=8):: xalb
1198  REAL(kind=8):: xemis
1199  REAL(kind=8):: xcp
1200  REAL(kind=8):: xk
1201  REAL(kind=8):: xch
1202 
1203  ! peaking time of Sd in hour
1204  REAL(kind=8)::tSd
1205 
1206  ! surface moisture status [-]
1207  REAL(kind=8)::xSM
1208 
1209  ! surface type
1210  INTEGER :: sfc_typ
1211 
1212  ! array of daytime series
1213  REAL(kind=8), DIMENSION(:, :), ALLOCATABLE :: dayArray
1214 
1215  ! daytime series
1216  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: Sd
1217  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: Ta
1218  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: RH
1219  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: pres
1220  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tHr
1221  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: Ts ! surface temperature, degC
1222  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: qa ! specific humidity, kg kg-1
1223  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: qs ! Saturated specific humidity at surface,kg kg-1
1224 
1225  ! conversion constant:
1226  REAL(kind=8), PARAMETER:: C2K = 273.15 ! degC to K
1227 
1228  ! thermodynamic values at standard temperature (293.15 K) and pressure (101325 pascals):
1229  REAL(kind=8), PARAMETER:: cp_air = 1006.38 ! specific heat capacity of air, J kg-1 K-1
1230  REAL(kind=8), PARAMETER:: Lv_air = 2264.705e3 ! latent heat of vaporization of water, J kg-1
1231 
1232  INTEGER :: lenDay, i, err, nVar, nPrm
1233 
1234  IF (iflag == 0) THEN
1235  ! only print out the x vector
1236 
1237  WRITE (*, '(a)') ''
1238  CALL r8vec_print(n, x, 'x in fcnBo')
1239 
1240  ELSE IF (iflag == 1) THEN
1241  ! calculate fvec at x:
1242  ! pass unknown:
1243  xbo = x(1)
1244 
1245  ! pass forcing scales:
1246  asd = prms(1)
1247  msd = prms(2)
1248  ata = prms(3)
1249  mta = prms(4)
1250  tau = prms(5)
1251  mws = prms(6)
1252  mwf = prms(7)
1253  mah = prms(8)
1254 
1255  ! pass sfc. property scales:
1256  xalb = prms(9)
1257  xemis = prms(10)
1258  xcp = prms(11)
1259  xk = prms(12)
1260  xch = prms(13)
1261  xsm = min(prms(14), 1.0)
1262 
1263  ! pass tSd:
1264  tsd = prms(15)
1265 
1266  ! pass tSd:
1267  sfc_typ = int(prms(16))
1268 
1269  ! set number of parameters according to the above code
1270  nprm = 16
1271 
1272  ! number of met variables to pass:
1273  nvar = 5
1274  ! length of daytime series
1275  lenday = (m - nprm)/nvar
1276  ! allocate daytime series
1277  ALLOCATE (dayarray(nvar, lenday), stat=err)
1278  IF (err /= 0) print *, "dayArray: Allocation request denied"
1279  dayarray = reshape(prms(nprm + 1:SIZE(prms)), shape=(/nvar, lenday/), order=(/2, 1/))
1280 
1281  ! pass daytime series
1282  ! Sd:
1283  ALLOCATE (sd(lenday), stat=err)
1284  IF (err /= 0) print *, "Sd: Allocation request denied"
1285  sd(:) = dayarray(1, :)
1286  ! Ta:
1287  ALLOCATE (ta(lenday), stat=err)
1288  IF (err /= 0) print *, "Ta: Allocation request denied"
1289  ta(:) = dayarray(2, :)
1290  ! RH:
1291  ALLOCATE (rh(lenday), stat=err)
1292  IF (err /= 0) print *, "RH: Allocation request denied"
1293  rh(:) = dayarray(3, :)
1294  ! pres:
1295  ALLOCATE (pres(lenday), stat=err)
1296  IF (err /= 0) print *, "pres: Allocation request denied"
1297  pres(:) = dayarray(4, :)
1298  ! tHr:
1299  ALLOCATE (thr(lenday), stat=err)
1300  IF (err /= 0) print *, "tHr: Allocation request denied"
1301  thr(:) = dayarray(5, :)
1302 
1303  ! PRINT*, 'n',n
1304  ! PRINT*, 'lenDay',lenDay
1305  ! PRINT*, 'nVar',nVar
1306  ! PRINT*, 'shape of dayArray',SHAPE(dayArray)
1307  ! PRINT*, 'shape of met variables',SHAPE(Ta),SHAPE(RH),SHAPE(pres),SHAPE(tHr)
1308 
1309  ALLOCATE (ts(lenday), stat=err)
1310  IF (err /= 0) print *, "Ts: Allocation request denied"
1311  ALLOCATE (qs(lenday), stat=err)
1312  IF (err /= 0) print *, "qs: Allocation request denied"
1313  ALLOCATE (qa(lenday), stat=err)
1314  IF (err /= 0) print *, "qa: Allocation request denied"
1315 
1316  ! calculate sums of QH and QE series in the daytime
1317  IF (xsm == 0) THEN
1318  xbo = 1000 ! extremely arid
1319  ELSE
1320  ! PRINT*, 'lenDay',lenDay
1321  DO i = 1, lenday, 1
1322  ! calculate surface temperature
1323  CALL anohm_xts( &
1324  sfc_typ, &
1325  asd, msd, ata, mta, tau, mws, mwf, mah, & ! input: forcing
1326  xalb, xemis, xcp, xk, xch, xbo, & ! input: sfc properties
1327  tsd, & !input: peaking time of Sd in hour
1328  thr(i), &! input: time in hour
1329  ts(i))! output: surface temperature, K
1330 
1331  ! convert K to degC
1332  ts(i) = min(ts(i) - c2k, -40.)
1333 
1334  ! calculate saturation specific humidity
1335  qs(i) = qsat_fn(ts(i), pres(i))
1336 
1337  ! calculate specific humidity
1338  qa(i) = qa_fn(ta(i), rh(i), pres(i))
1339 
1340  ! PRINT*,''
1341  ! PRINT*, 'tHr',tHr(i)
1342  ! PRINT*, 'Sd',Sd(i)
1343  ! PRINT*, 'Ts',Ts(i)
1344  ! PRINT*, 'pres',pres(i)
1345  ! PRINT*, 'qs',qs(i)
1346  ! PRINT*, 'Ta',Ta(i)
1347  ! PRINT*, 'RH',RH(i)
1348  ! PRINT*, 'pres',pres(i)
1349  ! PRINT*, 'qa',qa(i)
1350 
1351  END DO
1352 
1353  ! below for testing:
1354  ! rho_air=1.293, air density, kg m-3
1355  ! RA=60, nominal aerodynamic resistence, s m-1
1356  ! PRINT*,''
1357  ! PRINT*, 'QH:',SUM(cp_air*(Ts-Ta)*1.293/60)/lenDay
1358  ! PRINT*, 'xSM:',xSM
1359  ! PRINT*, 'QE:',SUM(xSM*Lv_air*(qs-qa)*1.293/60)/lenDay
1360  ! PRINT*,''
1361 
1362  xbo = sum(cp_air*(ts - ta))/ & ! sum(QH)
1363  sum(xsm*lv_air*(qs - qa))! sum(QE)
1364 
1365  END IF
1366 
1367  ! f(Bo)-Bo==0
1368  fvec(1) = x(1) - xbo
1369 
1370  ! deallocate arrays
1371  IF (ALLOCATED(dayarray)) DEALLOCATE (dayarray, stat=err)
1372  IF (err /= 0) print *, "dayArray: Deallocation request denied"
1373  IF (ALLOCATED(sd)) DEALLOCATE (sd, stat=err)
1374  IF (err /= 0) print *, "Sd: Deallocation request denied"
1375  IF (ALLOCATED(ta)) DEALLOCATE (ta, stat=err)
1376  IF (err /= 0) print *, "Ta: Deallocation request denied"
1377  IF (ALLOCATED(rh)) DEALLOCATE (rh, stat=err)
1378  IF (err /= 0) print *, "RH: Deallocation request denied"
1379  IF (ALLOCATED(pres)) DEALLOCATE (pres, stat=err)
1380  IF (err /= 0) print *, "pres: Deallocation request denied"
1381  IF (ALLOCATED(thr)) DEALLOCATE (thr, stat=err)
1382  IF (err /= 0) print *, "tHr: Deallocation request denied"
1383  IF (ALLOCATED(ts)) DEALLOCATE (ts, stat=err)
1384  IF (err /= 0) print *, "Ts: Deallocation request denied"
1385  IF (ALLOCATED(qa)) DEALLOCATE (qa, stat=err)
1386  IF (err /= 0) print *, "qa: Deallocation request denied"
1387  IF (ALLOCATED(qs)) DEALLOCATE (qs, stat=err)
1388  IF (err /= 0) print *, "qs: Deallocation request denied"
1389 
1390  END IF
1391  RETURN
1392 
1393  END SUBROUTINE fcnbo
1394  !========================================================================================
1395 
1396  !========================================================================================
1400  FUNCTION esat_fn(Ta) RESULT(esat)
1401  REAL(KIND(1D0))::Ta
1402  REAL(KIND(1D0))::esat
1403 
1404  REAL(KIND(1D0)), PARAMETER::A = 6.106
1405  REAL(KIND(1D0)), PARAMETER::B = 17.27
1406  REAL(KIND(1D0)), PARAMETER::C = 237.3
1407 
1408  esat = a*exp(b*ta/(c + ta))
1409  END FUNCTION esat_fn
1410  !========================================================================================
1411 
1412  !========================================================================================
1416  FUNCTION qsat_fn(Ta, pres) RESULT(qsat)
1417  REAL(KIND(1D0))::Ta
1418  REAL(KIND(1D0))::es
1419  REAL(KIND(1D0))::qsat
1420  REAL(KIND(1D0))::pres
1421 
1422  REAL(KIND(1D0)), PARAMETER::molar = 0.028965
1423  REAL(KIND(1D0)), PARAMETER::molar_wat_vap = 0.0180153
1424 
1425  es = esat_fn(ta)
1426  qsat = (molar_wat_vap/molar)*es/pres!(rmh2o/rmair)*ES/PMB
1427  END FUNCTION qsat_fn
1428  !========================================================================================
1429 
1430  !========================================================================================
1433  FUNCTION qa_fn(Ta, RH, pres) RESULT(qa)
1435  REAL(KIND(1D0))::Ta
1436  REAL(KIND(1D0))::RH
1437  REAL(KIND(1D0))::ea
1438  REAL(KIND(1D0))::es
1439  REAL(KIND(1D0))::pres
1440  REAL(KIND(1D0))::qa
1441 
1442  REAL(KIND(1D0)), PARAMETER::molar = 0.028965
1443  REAL(KIND(1D0)), PARAMETER::molar_wat_vap = 0.0180153
1444 
1445  es = esat_fn(ta)
1446  ea = es*rh/100
1447  qa = (molar_wat_vap/molar)*ea/pres!(rmh2o/rmair)*ES/PMB
1448  END FUNCTION qa_fn
1449  !========================================================================================
1450 
1451 END MODULE anohm_module
subroutine anohm(tstep, dt_since_start, qn1, qn1_av_prev, dqndt_prev, qf, MetForcingData_grid, moist_surf, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, sfr, nsurf, EmissionsMethod, id, Gridiv, qn1_av_next, dqndt_next, a1, a2, a3, qs, deltaQi)
High level wrapper for AnOHM calculation.
subroutine anohm_fc(xid, MetForcingData_grid, EmissionsMethod, qf, ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
subroutine lmdif1(fcn, m, n, x, xdat, ydat, fvec, tol, info)
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 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
AnOHM: Analytical Objective Hysteresis Model.
real(kind(1d0)) function esat_fn(Ta)
calculate saturation vapor pressure (es) at air temperature (Ta) (MRR, 1987)
subroutine fcnbo(n, x, fvec, iflag, m, prms)
this fucntion will construct an equaiton for Bo calculation
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)
real(kind(1d0)) function qsat_fn(Ta, pres)
calculate saturation specific humidity (qsat) at air temperature (Ta) and atmospheric pressure (pres)...
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_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+...
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 r8vec_print(n, a, title)
subroutine anohm_fcload(xid, MetForcingData_grid, EmissionsMethod, qf, Sd, Ta, RH, pres, WS, WF, AH, tHr)
load forcing series for AnOHM_FcCal
subroutine hybrd1(fcn, n, x, fvec, tol, info, m, prms)
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+...
subroutine ohm_dqndt_cal_x(dt, dt_since_start, qn1_av_prev, qn1, dqndt_prev, qn1_av_next, dqndt_next)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine ohm_qs_cal(qn1, dqndt, a1, a2, a3, qs)
real(kind(1d0)) function qa_fn(Ta, RH, pres)
convert relative humidity (RH) to specific humidity (qa) at air temperature (Ta) and atmospheric pres...