29 tstep, dt_since_start, &
30 qn1, qn_av_prev, dqndt_prev, qf, &
31 MetForcingData_grid, moist_surf, &
32 alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & ! input
33 sfr_surf, nsurf, EmissionsMethod, id, Gridiv, &
34 qn_av_next, dqndt_next, &
35 a1, a2, a3, qs, deltaQi)
38 INTEGER,
INTENT(in) :: tstep
39 INTEGER,
INTENT(in) :: dt_since_start
41 REAL(KIND(1D0)),
INTENT(in),
DIMENSION(:, :) :: MetForcingData_grid
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)
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
54 INTEGER,
INTENT(in) :: id
55 INTEGER,
INTENT(in) :: Gridiv
56 INTEGER,
INTENT(in) :: EmissionsMethod
57 INTEGER,
INTENT(in) :: nsurf
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
63 REAL(KIND(1D0)),
INTENT(out) :: dqndt_next
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)
75 INTEGER,
SAVE :: id_save
76 REAL(KIND(1D0)),
PARAMETER :: NotUsed = -55.5
77 INTEGER,
PARAMETER :: notUsedI = -55
82 REAL(KIND(1D0)),
DIMENSION(nsurf) :: xa1, xa2, xa3
94 idq = count(metforcingdata_grid(:, 2) == id .AND. &
95 metforcingdata_grid(:, 4) == 0 .AND. &
96 metforcingdata_grid(:, 15) > 0) &
111 CALL anohm_coef(is, xid, gridiv, metforcingdata_grid, moist_surf, emissionsmethod, qf, &
112 alb, emis, cpanohm, kkanohm, chanohm, &
113 xa1(is), xa2(is), xa3(is))
118 a1 = dot_product(xa1, sfr_surf)
119 a2 = dot_product(xa2, sfr_surf)
120 a3 = dot_product(xa3, sfr_surf)
128 CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn_av_prev, qn1, dqndt_prev, &
129 qn_av_next, dqndt_next)
132 CALL ohm_qs_cal(qn1, dqndt_prev, a1, a2, a3, qs)
135 CALL errorhint(21,
'SUEWS_AnOHM.f95: bad value for qn found during qs calculation.', qn1, notused, notusedi)
148 sfc_typ, xid, xgrid, & !input
149 MetForcingData_grid, moist, EmissionsMethod, qf, & !input
150 alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & ! input
156 INTEGER,
INTENT(in) :: sfc_typ
157 INTEGER,
INTENT(in) :: xid
158 INTEGER,
INTENT(in) :: xgrid
159 INTEGER,
INTENT(in) :: EmissionsMethod
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
171 REAL(KIND(1D0)),
INTENT(out) :: xa1
172 REAL(KIND(1D0)),
INTENT(out) :: xa2
173 REAL(KIND(1D0)),
INTENT(out) :: xa3
176 REAL(KIND(1D0)) :: ATs
177 REAL(KIND(1D0)) :: mTs
178 REAL(KIND(1D0)) :: gamma
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
187 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: Sd
188 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: Ta
189 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: RH
190 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: pres
191 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: WS
192 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: WF
193 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: AH
194 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: tHr
197 REAL(KIND(1D0)) :: xalb
198 REAL(KIND(1D0)) :: xemis
199 REAL(KIND(1D0)) :: xcp
200 REAL(KIND(1D0)) :: xk
201 REAL(KIND(1D0)) :: xch
202 REAL(KIND(1D0)) :: xBo
203 REAL(KIND(1D0)) :: xeta
204 REAL(KIND(1D0)) :: xmu
205 REAL(KIND(1D0)) :: xmoist
210 INTEGER,
SAVE :: id_save, grid_save
211 REAL(KIND(1D0)),
SAVE :: coeff_grid_day(7, 3) = -999.
217 IF (xid == id_save .AND. xgrid == grid_save)
THEN
220 xa1 = coeff_grid_day(sfc_typ, 1)
221 xa2 = coeff_grid_day(sfc_typ, 2)
222 xa3 = coeff_grid_day(sfc_typ, 3)
227 xid, metforcingdata_grid, emissionsmethod, qf, &
228 asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
232 xid, metforcingdata_grid, emissionsmethod, qf, &
233 sd, ta, rh, pres, ws, wf, ah, thr)
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)
247 sd, ta, rh, pres, thr, &
248 asd, msd, ata, mta, tau, mws, mwf, mah, &
249 xalb, xemis, xcp, xk, xch, xmoist, &
255 SELECT CASE (sfc_typ)
258 asd, msd, ata, mta, tau, mws, mwf, mah, &
259 xalb, xemis, xcp, xk, xch, xbo, &
260 xa1, xa2, xa3, ats, mts, gamma)
267 asd, msd, ata, mta, tau, mws, mwf, mah, &
268 xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, &
269 xa1, xa2, xa3, ats, mts, gamma)
277 coeff_grid_day(sfc_typ, :) = (/xa1, xa2, xa3/)
289 sfc_typ, & !input: surface type
290 ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
291 xalb, xemis, xcp, xk, xch, xBo, & ! input: sfc properties
292 tSd, & !input: peaking time of Sd in hour
293 xTHr, & ! input: time in hour
298 INTEGER,
INTENT(in) :: sfc_typ
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
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
321 REAL(KIND(1D0)),
INTENT(in) :: tSd
322 REAL(KIND(1D0)),
INTENT(in) :: xTHr
325 REAL(KIND(1D0)),
INTENT(out) :: xTs
333 REAL(KIND(1D0)),
PARAMETER :: PI = atan(1.0)*4
334 REAL(KIND(1D0)),
PARAMETER :: OMEGA = 2*pi/(24*60*60)
336 SELECT CASE (sfc_typ)
339 asd, msd, ata, mta, tau, mws, mwf, mah, &
340 xalb, xemis, xcp, xk, xch, xbo, &
341 xa1, xa2, xa3, ats, mts, gamma)
348 asd, msd, ata, mta, tau, mws, mwf, mah, &
349 xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, &
350 xa1, xa2, xa3, ats, mts, gamma)
355 xts = ats*sin(omega*(xthr - tsd + 6)*3600 - gamma) + mts
362 ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
363 xalb, xemis, xcp, xk, xch, xBo, & ! input: sfc properties
364 xa1, xa2, xa3, ATs, mTs, gamma)
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
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
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
394 REAL(KIND(1D0)),
PARAMETER :: SIGMA = 5.67e-8
395 REAL(KIND(1D0)),
PARAMETER :: PI = atan(1.0)*4
396 REAL(KIND(1D0)),
PARAMETER :: OMEGA = 2*pi/(24*60*60)
399 REAL(KIND(1D0)) :: beta
400 REAL(KIND(1D0)) :: f, fL, fT
401 REAL(KIND(1D0)) :: lambda
402 REAL(KIND(1D0)) :: delta, m, n
403 REAL(KIND(1D0)) :: ms, ns
404 REAL(KIND(1D0)) :: ceta, cphi
405 REAL(KIND(1D0)) :: eta, phi, xlag
406 REAL(KIND(1D0)) :: xx1, xx2, xx3, xchWS
436 IF (abs(xbo) < 0.001)
THEN
443 f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
445 fl = 4*sigma*xemis*mta**3
447 ft = (1 + beta)*xchws
451 delta = sqrt(.5*(mwf**2 + sqrt(mwf**4 + 16*lambda**2*omega**2)))
453 m = (2*lambda)/(delta + mwf)
458 mts = (msd*(1 - xalb)/f) + mta
464 xx2 = (1 - xalb)*asd + f*cos(tau)*ata
466 gamma = atan(ns/ms) + atan(xx1/xx2)
468 ats = -(sin(tau)*ata)/(ns*cos(gamma) - ms*sin(gamma))
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
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)
483 xx1 = m*cos(gamma) - n*sin(gamma)
484 xx2 = m*sin(gamma) + n*cos(gamma)
487 xx1 = xk**2*(m**2 + n**2)*ats**2
496 xa1 = (ceta*cos(xlag))/cphi
499 xa2 = (ceta*sin(xlag))/(omega*cphi)
503 xa3 = -xa1*(ft/f)*(msd*(1 - xalb)) - mah
519 ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
520 xalb, xemis, xcp, xk, xch, xBo, xeta, xmu, & ! input: sfc properties
521 xa1, xa2, xa3, ATs, mTs, gamma)
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
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
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
554 REAL(KIND(1D0)),
PARAMETER :: SIGMA = 5.67e-8
555 REAL(KIND(1D0)),
PARAMETER :: PI = atan(1.0)*4
556 REAL(KIND(1D0)),
PARAMETER :: OMEGA = 2*pi/(24*60*60)
559 REAL(KIND(1D0)) :: beta
560 REAL(KIND(1D0)) :: f, fL, fT
561 REAL(KIND(1D0)) :: lambda, calb
562 REAL(KIND(1D0)) :: delta
564 REAL(KIND(1D0)) :: xm, xn
566 REAL(KIND(1D0)) :: phi
568 REAL(KIND(1D0)) :: czeta, ctheta
569 REAL(KIND(1D0)) :: zeta, theta, xlag
570 REAL(KIND(1D0)) :: xx1, xx2, xx3
571 REAL(KIND(1D0)) :: kappa
572 REAL(KIND(1D0)) :: dtau, dpsi, dphi
573 REAL(KIND(1D0)) :: cdtau, cdpsi, cdphi
574 REAL(KIND(1D0)) :: xxT, xxkappa, xxdltphi, xchWS
578 REAL(KIND(1D0)) :: dummy
586 kappa = sqrt(xcp*omega/(2*xk))
591 f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
592 fl = 4*sigma*xemis*mta**3
593 ft = (1 + beta)*xchws
597 lambda = sqrt(xm**2 + xn**2)
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)
608 mts = (msd*(1 - xalb + xeta)/f) + mta
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)
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)
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)
628 ctheta = sqrt(xx1**2 + xx2**2)
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)
636 xx1 = xxt*sin(pi/4 - delta) + xxkappa*sin(phi + dpsi)
637 xx2 = xxt*sin(pi/4 + delta) - xxkappa*sin(phi - dpsi)
640 xx1 = 2*sqrt(2.)*xxkappa*xxt*xxdltphi
641 xx2 = (1 - cos(2*dpsi)*cos(2*phi))*xxkappa**2
643 czeta = xk/lambda*sqrt(xx1 + xx2 + xx3)
648 xa1 = (czeta*cos(xlag))/ctheta
652 xa2 = (czeta*sin(xlag))/(omega*ctheta)
656 xa3 = msd*(xalb - 1)*(xeta + (ft - fl*xeta)/f*xa1)
666 xid, MetForcingData_grid, EmissionsMethod, qf, & ! input
667 ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
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
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
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
695 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: WF
696 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: AH
697 REAL(KIND(1D0)),
DIMENSION(:),
ALLOCATABLE :: tHr
700 CALL anohm_fcload(xid, metforcingdata_grid, emissionsmethod, qf, sd, ta, rh, pres, ws, wf, ah, thr)
702 CALL anohm_fccal(sd, ta, ws, wf, ah, thr, asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
716 xid, MetForcingData_grid, EmissionsMethod, qf, & ! input
717 Sd, Ta, RH, pres, WS, WF, AH, tHr)
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
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
738 REAL(KIND(1D0)),
DIMENSION(:, :),
ALLOCATABLE :: subMet
741 INTEGER :: lenMetData, nVar
743 LOGICAL,
ALLOCATABLE :: metMask(:)
746 IF (
ALLOCATED(metmask))
DEALLOCATE (metmask, stat=err)
747 ALLOCATE (metmask(
SIZE(metforcingdata_grid, dim=1)))
748 metmask = (metforcingdata_grid(:, 2) == xid &
749 .AND. metforcingdata_grid(:, 4) == 0)
752 lenmetdata = count(metmask)
758 IF (
ALLOCATED(submet))
DEALLOCATE (submet, stat=err)
759 ALLOCATE (submet(lenmetdata, nvar))
760 submet = reshape(pack(metforcingdata_grid(:, (/3, &
761 15, 12, 11, 13, 10, 12, 9/)), &
763 spread(metmask, dim=2, ncopies=nvar)), &
764 (/lenmetdata, nvar/))
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))
792 IF (emissionsmethod == 0)
THEN
807 Sd, Ta, WS, WF, AH, tHr, & ! input
808 ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
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
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
832 REAL(KIND(1D0)),
PARAMETER :: PI = atan(1.0)*4
833 REAL(KIND(1D0)),
PARAMETER :: C2K = 273.15
836 REAL(KIND(1D0)),
ALLOCATABLE :: tHrDay(:)
837 REAL(KIND(1D0)),
ALLOCATABLE :: selX(:)
840 INTEGER :: err, lenDay
841 LOGICAL,
DIMENSION(:),
ALLOCATABLE :: SdMask
843 ALLOCATE (sdmask(
SIZE(sd, dim=1)), stat=err)
844 IF (err /= 0) print *,
"SdMask: Allocation request denied"
846 lenday = count(sdmask)
850 ALLOCATE (thrday(lenday), stat=err)
851 IF (err /= 0) print *,
"tHrDay: Allocation request denied"
852 thrday = pack(thr, mask=sdmask)
854 ALLOCATE (selx(lenday), stat=err)
855 IF (err /= 0) print *,
"selX: Allocation request denied"
859 selx = pack(sd, mask=sdmask)
860 asd = (maxval(selx) - minval(selx))/2
861 msd = sum(selx)/lenday
878 IF (msd + asd < maxval(selx))
THEN
882 msd = maxval(selx) - asd
889 selx = pack(ta, mask=sdmask)
890 ata = (maxval(selx) - minval(selx))/2
891 mta = sum(selx)/lenday
895 IF (mta < 60) mta = mta + c2k
908 tau = (tta - tsd)/24*2*pi
912 selx = pack(ws, mask=sdmask)
913 mws = sum(selx)/lenday
915 selx = pack(wf, mask=sdmask)
916 mwf = sum(selx)/lenday
918 selx = pack(ah, mask=sdmask)
919 mah = sum(selx)/lenday
923 IF (
ALLOCATED(sdmask))
DEALLOCATE (sdmask, stat=err)
924 IF (err /= 0) print *,
"SdMask: Deallocation request denied"
926 IF (
ALLOCATED(thrday))
DEALLOCATE (thrday, stat=err)
927 IF (err /= 0) print *,
"tHrDay: Deallocation request denied"
929 IF (
ALLOCATED(selx))
DEALLOCATE (selx, stat=err)
930 IF (err /= 0) print *,
"selX: Deallocation request denied"
944 REAL(KIND(1D0)),
INTENT(in) :: tHr(:)
945 REAL(KIND(1D0)),
INTENT(in) :: obs(:)
948 REAL(KIND(1D0)),
INTENT(out) :: amp
949 REAL(KIND(1D0)),
INTENT(out) :: mean
950 REAL(KIND(1D0)),
INTENT(out) :: tpeak
952 INTEGER :: m, n, info, err
955 REAL(KIND(1D0)),
ALLOCATABLE :: fvec(:), x(:)
957 REAL(KIND(1D0)) :: tol = 0.00001d+00
963 ALLOCATE (fvec(m), stat=err)
964 IF (err /= 0) print *,
"fvec: Allocation request denied"
966 ALLOCATE (x(n), stat=err)
967 IF (err /= 0) print *,
"x: Allocation request denied"
970 x = (/mean, amp, tpeak/)
973 CALL lmdif1(
fsin, m, n, x, thr, obs, fvec, tol, info)
989 tpeak = x(3) - 12 + 6 + 24
990 tpeak = mod(tpeak, 24.)
994 tpeak = mod(tpeak, 24.)
1003 SUBROUTINE fsin(m, n, x, xdat, ydat, fvec, iflag)
1010 REAL(kind=8) fvec(m), &
1013 INTEGER(kind=4) iflag, i
1015 REAL(KIND(1D0)),
PARAMETER :: PI = atan(1.0)*4
1017 IF (iflag == 0)
THEN
1021 WRITE (*,
'(g14.6)') x(i)
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)
1036 sfc_typ, & ! surface type
1037 Sd, Ta, RH, pres, tHr, & ! input: forcing
1038 ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing
1039 xalb, xemis, xcp, xk, xch, xSM, & ! input: sfc properties
1040 tSd, & ! input: peaking time of Sd in hour
1046 INTEGER,
INTENT(in) :: sfc_typ
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
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
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
1075 REAL(kind=8), intent(out) :: xbo
1079 REAL(kind=8), allocatable :: x(:), fvec(:), prms(:)
1081 INTEGER :: lenDay, n, m, info, err, nVar, nPrm
1082 LOGICAL,
DIMENSION(:),
ALLOCATABLE :: maskDay
1083 REAL(kind=8) :: tol = 1e-20
1088 ALLOCATE (maskday(
SIZE(sd)), stat=err)
1091 lenday =
SIZE(pack(sd, mask=maskday), dim=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"
1114 m = nprm + nvar*lenday
1116 ALLOCATE (prms(m), stat=err)
1117 IF (err /= 0) print *,
"prms: Allocation request denied"
1141 prms(16) = sfc_typ*1.0
1144 prms(nprm + 1:m) = pack((/sd, ta, rh, pres, thr/), &
1145 mask=pack(spread(maskday, dim=2, ncopies=nvar), .true.))
1151 CALL hybrd1(
fcnbo, n, x, fvec, tol, info, m, prms)
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"
1168 SUBROUTINE fcnbo(n, x, fvec, iflag, m, prms)
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)
1179 REAL(kind=8) :: prms(m)
1195 REAL(kind=8) :: xalb
1196 REAL(kind=8) :: xemis
1211 REAL(kind=8), dimension(:, :),
ALLOCATABLE :: dayarray
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
1220 REAL(kind=8), dimension(:),
ALLOCATABLE :: qa
1221 REAL(kind=8), dimension(:),
ALLOCATABLE :: qs
1224 REAL(kind=8), parameter :: c2k = 273.15
1227 REAL(kind=8), parameter :: cp_air = 1006.38
1228 REAL(kind=8), parameter :: lv_air = 2264.705e3
1230 INTEGER :: lenDay, i, err, nVar, nPrm
1232 IF (iflag == 0)
THEN
1238 ELSE IF (iflag == 1)
THEN
1259 xsm = min(prms(14), 1.0)
1265 sfc_typ = int(prms(16))
1273 lenday = (m - nprm)/nvar
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/))
1281 ALLOCATE (sd(lenday), stat=err)
1282 IF (err /= 0) print *,
"Sd: Allocation request denied"
1283 sd(:) = dayarray(1, :)
1285 ALLOCATE (ta(lenday), stat=err)
1286 IF (err /= 0) print *,
"Ta: Allocation request denied"
1287 ta(:) = dayarray(2, :)
1289 ALLOCATE (rh(lenday), stat=err)
1290 IF (err /= 0) print *,
"RH: Allocation request denied"
1291 rh(:) = dayarray(3, :)
1293 ALLOCATE (pres(lenday), stat=err)
1294 IF (err /= 0) print *,
"pres: Allocation request denied"
1295 pres(:) = dayarray(4, :)
1297 ALLOCATE (thr(lenday), stat=err)
1298 IF (err /= 0) print *,
"tHr: Allocation request denied"
1299 thr(:) = dayarray(5, :)
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"
1323 asd, msd, ata, mta, tau, mws, mwf, mah, &
1324 xalb, xemis, xcp, xk, xch, xbo, &
1330 ts(i) = min(ts(i) - c2k, -40.)
1333 qs(i) =
qsat_fn(ts(i), pres(i))
1336 qa(i) =
qa_fn(ta(i), rh(i), pres(i))
1360 xbo = sum(cp_air*(ts - ta))/ &
1361 sum(xsm*lv_air*(qs - qa))
1366 fvec(1) = x(1) - xbo
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"
1391 END SUBROUTINE fcnbo
1399 REAL(kind(1d0)) :: ta
1400 REAL(kind(1d0)) :: esat
1402 REAL(kind(1d0)),
PARAMETER :: a = 6.106
1403 REAL(kind(1d0)),
PARAMETER :: b = 17.27
1404 REAL(kind(1d0)),
PARAMETER :: c = 237.3
1406 esat = a*exp(b*ta/(c + ta))
1415 REAL(kind(1d0)) :: ta
1416 REAL(kind(1d0)) :: es
1417 REAL(kind(1d0)) :: qsat
1418 REAL(kind(1d0)) :: pres
1420 REAL(kind(1d0)),
PARAMETER :: molar = 0.028965
1421 REAL(kind(1d0)),
PARAMETER :: molar_wat_vap = 0.0180153
1424 qsat = (molar_wat_vap/molar)*es/pres
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
1440 REAL(kind(1d0)),
PARAMETER :: molar = 0.028965
1441 REAL(kind(1d0)),
PARAMETER :: molar_wat_vap = 0.0180153
1445 qa = (molar_wat_vap/molar)*ea/pres
AnOHM: Analytical Objective Hysteresis Model.
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 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_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 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_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 fcnbo(n, x, fvec, iflag, m, prms)
this fucntion will construct an equaiton for Bo calculation
real(kind(1d0)) function qa_fn(ta, rh, pres)
convert relative humidity (RH) to specific humidity (qa) at air temperature (Ta) and atmospheric pres...
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.
real(kind(1d0)) function qsat_fn(ta, pres)
calculate saturation specific humidity (qsat) at air temperature (Ta) and atmospheric pressure (pres)...
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,...
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)
real(kind(1d0)) function esat_fn(ta)
calculate saturation vapor pressure (es) at air temperature (Ta) (MRR, 1987)
subroutine anohm_fcload(xid, metforcingdata_grid, emissionsmethod, qf, sd, ta, rh, pres, ws, wf, ah, thr)
load forcing series for AnOHM_FcCal
subroutine anohm_fc(xid, metforcingdata_grid, emissionsmethod, qf, asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
subroutine errorhint(errh, problemfile, value, value2, valuei)
subroutine ohm_qs_cal(qn1, dqndt, a1, a2, a3, qs)
subroutine ohm_dqndt_cal_x(dt, dt_since_start, qn1_av_prev, qn1, dqndt_prev, qn1_av_next, dqndt_next)
subroutine lmdif1(fcn, m, n, x, xdat, ydat, fvec, tol, info)
subroutine hybrd1(fcn, n, x, fvec, tol, info, m, prms)
subroutine r8vec_print(n, a, title)