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)
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(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)::qn1_av_prev
61 REAL(KIND(1d0)),
INTENT(out)::qn1_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
96 idq = count(metforcingdata_grid(:, 2) == id .AND. &
97 metforcingdata_grid(:, 4) == 0 .AND. &
98 metforcingdata_grid(:, 15) > 0) &
113 CALL anohm_coef(is, xid, gridiv, metforcingdata_grid, moist_surf, emissionsmethod, qf, &
114 alb, emis, cpanohm, kkanohm, chanohm, &
115 xa1(is), xa2(is), xa3(is))
120 a1 = dot_product(xa1, sfr)
121 a2 = dot_product(xa2, sfr)
122 a3 = dot_product(xa3, sfr)
130 CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn1_av_prev, qn1, dqndt_prev, &
131 qn1_av_next, dqndt_next)
134 CALL ohm_qs_cal(qn1, dqndt_prev, a1, a2, a3, qs)
137 CALL errorhint(21,
'SUEWS_AnOHM.f95: bad value for qn found during qs calculation.', qn1, notused, notusedi)
150 sfc_typ, xid, xgrid, &!input
151 MetForcingData_grid, moist, EmissionsMethod, qf, & !input
152 alb, emis, cpAnOHM, kkAnOHM, chAnOHM, &! input
158 INTEGER,
INTENT(in):: sfc_typ
159 INTEGER,
INTENT(in):: xid
160 INTEGER,
INTENT(in):: xgrid
161 INTEGER,
INTENT(in):: EmissionsMethod
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
173 REAL(KIND(1d0)),
INTENT(out) :: xa1
174 REAL(KIND(1d0)),
INTENT(out) :: xa2
175 REAL(KIND(1d0)),
INTENT(out) :: xa3
178 REAL(KIND(1d0)):: ATs
179 REAL(KIND(1d0)):: mTs
180 REAL(KIND(1d0)):: gamma
183 REAL(KIND(1d0))::ASd, mSd, tSd
184 REAL(KIND(1d0))::ATa, mTa, tTa
186 REAL(KIND(1d0))::mWS, mWF, mAH
189 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::Sd
190 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::Ta
191 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::RH
192 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::pres
193 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::WS
194 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::WF
195 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::AH
196 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE::tHr
199 REAL(KIND(1d0)) ::xalb
200 REAL(KIND(1d0)) ::xemis
201 REAL(KIND(1d0)) ::xcp
203 REAL(KIND(1d0)) ::xch
204 REAL(KIND(1d0)) ::xBo
205 REAL(KIND(1d0)) ::xeta
206 REAL(KIND(1d0)) ::xmu
207 REAL(KIND(1d0)) ::xmoist
212 INTEGER,
SAVE :: id_save, grid_save
213 REAL(KIND(1d0)),
SAVE:: coeff_grid_day(7, 3) = -999.
219 IF (xid == id_save .AND. xgrid == grid_save)
THEN 222 xa1 = coeff_grid_day(sfc_typ, 1)
223 xa2 = coeff_grid_day(sfc_typ, 2)
224 xa3 = coeff_grid_day(sfc_typ, 3)
229 xid, metforcingdata_grid, emissionsmethod, qf, &
230 asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
234 xid, metforcingdata_grid, emissionsmethod, qf, &
235 sd, ta, rh, pres, ws, wf, ah, thr)
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)
249 sd, ta, rh, pres, thr, &
250 asd, msd, ata, mta, tau, mws, mwf, mah, &
251 xalb, xemis, xcp, xk, xch, xmoist, &
257 SELECT CASE (sfc_typ)
260 asd, msd, ata, mta, tau, mws, mwf, mah, &
261 xalb, xemis, xcp, xk, xch, xbo, &
262 xa1, xa2, xa3, ats, mts, gamma)
269 asd, msd, ata, mta, tau, mws, mwf, mah, &
270 xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, &
271 xa1, xa2, xa3, ats, mts, gamma)
279 coeff_grid_day(sfc_typ, :) = (/xa1, xa2, xa3/)
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
300 INTEGER,
INTENT(in):: sfc_typ
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
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
323 REAL(KIND(1d0)),
INTENT(in):: tSd
324 REAL(KIND(1d0)),
INTENT(in):: xTHr
327 REAL(KIND(1d0)),
INTENT(out) :: xTs
335 REAL(KIND(1d0)),
PARAMETER :: PI = atan(1.0)*4
336 REAL(KIND(1d0)),
PARAMETER :: OMEGA = 2*pi/(24*60*60)
338 SELECT CASE (sfc_typ)
341 asd, msd, ata, mta, tau, mws, mwf, mah, &
342 xalb, xemis, xcp, xk, xch, xbo, &
343 xa1, xa2, xa3, ats, mts, gamma)
350 asd, msd, ata, mta, tau, mws, mwf, mah, &
351 xalb, xemis, xcp, xk, xch, xbo, xeta, xmu, &
352 xa1, xa2, xa3, ats, mts, gamma)
357 xts = ats*sin(omega*(xthr - tsd + 6)*3600 - gamma) + mts
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)
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
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
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
396 REAL(KIND(1d0)),
PARAMETER :: SIGMA = 5.67e-8
397 REAL(KIND(1d0)),
PARAMETER :: PI = atan(1.0)*4
398 REAL(KIND(1d0)),
PARAMETER :: OMEGA = 2*pi/(24*60*60)
401 REAL(KIND(1d0)) :: beta
402 REAL(KIND(1d0)) :: f, fL, fT
403 REAL(KIND(1d0)) :: lambda
404 REAL(KIND(1d0)) :: delta, m, n
405 REAL(KIND(1d0)) :: ms, ns
406 REAL(KIND(1d0)) :: ceta, cphi
407 REAL(KIND(1d0)) :: eta, phi, xlag
408 REAL(KIND(1d0)) :: xx1, xx2, xx3, xchWS
438 IF (abs(xbo) < 0.001)
THEN 445 f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
447 fl = 4*sigma*xemis*mta**3
449 ft = (1 + beta)*xchws
453 delta = sqrt(.5*(mwf**2 + sqrt(mwf**4 + 16*lambda**2*omega**2)))
455 m = (2*lambda)/(delta + mwf)
460 mts = (msd*(1 - xalb)/f) + mta
466 xx2 = (1 - xalb)*asd + f*cos(tau)*ata
468 gamma = atan(ns/ms) + atan(xx1/xx2)
470 ats = -(sin(tau)*ata)/(ns*cos(gamma) - ms*sin(gamma))
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
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)
485 xx1 = m*cos(gamma) - n*sin(gamma)
486 xx2 = m*sin(gamma) + n*cos(gamma)
489 xx1 = xk**2*(m**2 + n**2)*ats**2
498 xa1 = (ceta*cos(xlag))/cphi
501 xa2 = (ceta*sin(xlag))/(omega*cphi)
505 xa3 = -xa1*(ft/f)*(msd*(1 - xalb)) - mah
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)
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
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
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
556 REAL(KIND(1d0)),
PARAMETER :: SIGMA = 5.67e-8
557 REAL(KIND(1d0)),
PARAMETER :: PI = atan(1.0)*4
558 REAL(KIND(1d0)),
PARAMETER :: OMEGA = 2*pi/(24*60*60)
561 REAL(KIND(1d0)) :: beta
562 REAL(KIND(1d0)) :: f, fL, fT
563 REAL(KIND(1d0)) :: lambda, calb
564 REAL(KIND(1d0)) :: delta
566 REAL(KIND(1d0)) :: xm, xn
568 REAL(KIND(1d0)) :: phi
570 REAL(KIND(1d0)) :: czeta, ctheta
571 REAL(KIND(1d0)) :: zeta, theta, xlag
572 REAL(KIND(1d0)) :: xx1, xx2, xx3
573 REAL(KIND(1d0)) :: kappa
574 REAL(KIND(1d0)) :: dtau, dpsi, dphi
575 REAL(KIND(1d0)) :: cdtau, cdpsi, cdphi
576 REAL(KIND(1d0)) :: xxT, xxkappa, xxdltphi, xchWS
580 REAL(KIND(1d0)) :: dummy
588 kappa = sqrt(xcp*omega/(2*xk))
593 f = ((1 + beta)*xchws + 4*sigma*xemis*mta**3)
594 fl = 4*sigma*xemis*mta**3
595 ft = (1 + beta)*xchws
599 lambda = sqrt(xm**2 + xn**2)
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)
610 mts = (msd*(1 - xalb + xeta)/f) + mta
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)
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)
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)
630 ctheta = sqrt(xx1**2 + xx2**2)
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)
638 xx1 = xxt*sin(pi/4 - delta) + xxkappa*sin(phi + dpsi)
639 xx2 = xxt*sin(pi/4 + delta) - xxkappa*sin(phi - dpsi)
642 xx1 = 2*sqrt(2.)*xxkappa*xxt*xxdltphi
643 xx2 = (1 - cos(2*dpsi)*cos(2*phi))*xxkappa**2
645 czeta = xk/lambda*sqrt(xx1 + xx2 + xx3)
650 xa1 = (czeta*cos(xlag))/ctheta
654 xa2 = (czeta*sin(xlag))/(omega*ctheta)
658 xa3 = msd*(xalb - 1)*(xeta + (ft - fl*xeta)/f*xa1)
668 xid, MetForcingData_grid, EmissionsMethod, qf, & ! input
669 ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
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
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
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
697 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE:: WF
698 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE:: AH
699 REAL(KIND(1d0)),
DIMENSION(:),
ALLOCATABLE:: tHr
702 CALL anohm_fcload(xid, metforcingdata_grid, emissionsmethod, qf, sd, ta, rh, pres, ws, wf, ah, thr)
704 CALL anohm_fccal(sd, ta, ws, wf, ah, thr, asd, msd, tsd, ata, mta, tta, tau, mws, mwf, mah)
718 xid, MetForcingData_grid, EmissionsMethod, qf, & ! input
719 Sd, Ta, RH, pres, WS, WF, AH, tHr)
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
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
740 REAL(KIND(1d0)),
DIMENSION(:, :),
ALLOCATABLE :: subMet
743 INTEGER :: lenMetData, nVar
745 LOGICAL,
ALLOCATABLE :: metMask(:)
748 IF (
ALLOCATED(metmask))
DEALLOCATE (metmask, stat=err)
749 ALLOCATE (metmask(
SIZE(metforcingdata_grid, dim=1)))
750 metmask = (metforcingdata_grid(:, 2) == xid &
751 .AND. metforcingdata_grid(:, 4) == 0)
754 lenmetdata = count(metmask)
760 IF (
ALLOCATED(submet))
DEALLOCATE (submet, stat=err)
761 ALLOCATE (submet(lenmetdata, nvar))
762 submet = reshape(pack(metforcingdata_grid(:, (/3, &
763 15, 12, 11, 13, 10, 12, 9/)), &
765 spread(metmask, dim=2, ncopies=nvar)), &
766 (/lenmetdata, nvar/))
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))
794 IF (emissionsmethod == 0)
THEN 809 Sd, Ta, WS, WF, AH, tHr, & ! input
810 ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH)
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
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
834 REAL(KIND(1d0)),
PARAMETER :: PI = atan(1.0)*4
835 REAL(KIND(1d0)),
PARAMETER :: C2K = 273.15
838 REAL(KIND(1d0)),
ALLOCATABLE :: tHrDay(:)
839 REAL(KIND(1d0)),
ALLOCATABLE :: selX(:)
842 INTEGER :: err, lenDay
843 LOGICAL,
DIMENSION(:),
ALLOCATABLE::SdMask
845 ALLOCATE (sdmask(
SIZE(sd, dim=1)), stat=err)
846 IF (err /= 0) print *,
"SdMask: Allocation request denied" 848 lenday = count(sdmask)
852 ALLOCATE (thrday(lenday), stat=err)
853 IF (err /= 0) print *,
"tHrDay: Allocation request denied" 854 thrday = pack(thr, mask=sdmask)
856 ALLOCATE (selx(lenday), stat=err)
857 IF (err /= 0) print *,
"selX: Allocation request denied" 861 selx = pack(sd, mask=sdmask)
862 asd = (maxval(selx) - minval(selx))/2
863 msd = sum(selx)/lenday
880 IF (msd + asd < maxval(selx))
THEN 884 msd = maxval(selx) - asd
891 selx = pack(ta, mask=sdmask)
892 ata = (maxval(selx) - minval(selx))/2
893 mta = sum(selx)/lenday
897 IF (mta < 60) mta = mta + c2k
910 tau = (tta - tsd)/24*2*pi
914 selx = pack(ws, mask=sdmask)
915 mws = sum(selx)/lenday
917 selx = pack(wf, mask=sdmask)
918 mwf = sum(selx)/lenday
920 selx = pack(ah, mask=sdmask)
921 mah = sum(selx)/lenday
925 IF (
ALLOCATED(sdmask))
DEALLOCATE (sdmask, stat=err)
926 IF (err /= 0) print *,
"SdMask: Deallocation request denied" 928 IF (
ALLOCATED(thrday))
DEALLOCATE (thrday, stat=err)
929 IF (err /= 0) print *,
"tHrDay: Deallocation request denied" 931 IF (
ALLOCATED(selx))
DEALLOCATE (selx, stat=err)
932 IF (err /= 0) print *,
"selX: Deallocation request denied" 946 REAL(KIND(1d0)),
INTENT(in) :: tHr(:)
947 REAL(KIND(1d0)),
INTENT(in) :: obs(:)
950 REAL(KIND(1d0)),
INTENT(out) :: amp
951 REAL(KIND(1d0)),
INTENT(out) :: mean
952 REAL(KIND(1d0)),
INTENT(out) :: tpeak
954 INTEGER :: m, n, info, err
957 REAL(KIND(1d0)),
ALLOCATABLE:: fvec(:), x(:)
959 REAL(KIND(1d0)):: tol = 0.00001d+00
965 ALLOCATE (fvec(m), stat=err)
966 IF (err /= 0) print *,
"fvec: Allocation request denied" 968 ALLOCATE (x(n), stat=err)
969 IF (err /= 0) print *,
"x: Allocation request denied" 972 x = (/mean, amp, tpeak/)
975 CALL lmdif1(
fsin, m, n, x, thr, obs, fvec, tol, info)
991 tpeak = x(3) - 12 + 6 + 24
992 tpeak = mod(tpeak, 24.)
996 tpeak = mod(tpeak, 24.)
1005 SUBROUTINE fsin(m, n, x, xdat, ydat, fvec, iflag)
1012 REAL(kind=8) fvec(m), &
1015 INTEGER(kind=4) iflag, i
1017 REAL(KIND(1d0)),
PARAMETER :: PI = atan(1.0)*4
1019 IF (iflag == 0)
THEN 1023 WRITE (*,
'(g14.6)') x(i)
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)
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
1048 INTEGER,
INTENT(in) :: sfc_typ
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
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
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
1077 REAL(kind=8),
INTENT(out) :: xBo
1081 REAL(kind=8),
ALLOCATABLE :: x(:), fvec(:), prms(:)
1083 INTEGER :: lenDay, n, m, info, err, nVar, nPrm
1084 LOGICAL,
DIMENSION(:),
ALLOCATABLE :: maskDay
1085 REAL(kind=8) :: tol = 1e-20
1090 ALLOCATE (maskday(
SIZE(sd)), stat=err)
1093 lenday =
SIZE(pack(sd, mask=maskday), dim=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" 1116 m = nprm + nvar*lenday
1118 ALLOCATE (prms(m), stat=err)
1119 IF (err /= 0) print *,
"prms: Allocation request denied" 1143 prms(16) = sfc_typ*1.0
1146 prms(nprm + 1:m) = pack((/sd, ta, rh, pres, thr/), &
1147 mask=pack(spread(maskday, dim=2, ncopies=nvar), .true.))
1153 CALL hybrd1(
fcnbo, n, x, fvec, tol, info, m, prms)
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" 1170 SUBROUTINE fcnbo(n, x, fvec, iflag, m, prms)
1178 INTEGER(kind=4):: iflag
1179 REAL(kind=8):: fvec(n)
1181 REAL(kind=8):: prms(m)
1198 REAL(kind=8):: xemis
1213 REAL(kind=8),
DIMENSION(:, :),
ALLOCATABLE :: dayArray
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
1222 REAL(kind=8),
DIMENSION(:),
ALLOCATABLE :: qa
1223 REAL(kind=8),
DIMENSION(:),
ALLOCATABLE :: qs
1226 REAL(kind=8),
PARAMETER:: C2K = 273.15
1229 REAL(kind=8),
PARAMETER:: cp_air = 1006.38
1230 REAL(kind=8),
PARAMETER:: Lv_air = 2264.705e3
1232 INTEGER :: lenDay, i, err, nVar, nPrm
1234 IF (iflag == 0)
THEN 1240 ELSE IF (iflag == 1)
THEN 1261 xsm = min(prms(14), 1.0)
1267 sfc_typ = int(prms(16))
1275 lenday = (m - nprm)/nvar
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/))
1283 ALLOCATE (sd(lenday), stat=err)
1284 IF (err /= 0) print *,
"Sd: Allocation request denied" 1285 sd(:) = dayarray(1, :)
1287 ALLOCATE (ta(lenday), stat=err)
1288 IF (err /= 0) print *,
"Ta: Allocation request denied" 1289 ta(:) = dayarray(2, :)
1291 ALLOCATE (rh(lenday), stat=err)
1292 IF (err /= 0) print *,
"RH: Allocation request denied" 1293 rh(:) = dayarray(3, :)
1295 ALLOCATE (pres(lenday), stat=err)
1296 IF (err /= 0) print *,
"pres: Allocation request denied" 1297 pres(:) = dayarray(4, :)
1299 ALLOCATE (thr(lenday), stat=err)
1300 IF (err /= 0) print *,
"tHr: Allocation request denied" 1301 thr(:) = dayarray(5, :)
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" 1325 asd, msd, ata, mta, tau, mws, mwf, mah, &
1326 xalb, xemis, xcp, xk, xch, xbo, &
1332 ts(i) = min(ts(i) - c2k, -40.)
1335 qs(i) =
qsat_fn(ts(i), pres(i))
1338 qa(i) =
qa_fn(ta(i), rh(i), pres(i))
1362 xbo = sum(cp_air*(ts - ta))/ &
1363 sum(xsm*lv_air*(qs - qa))
1368 fvec(1) = x(1) - xbo
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" 1393 END SUBROUTINE fcnbo 1400 FUNCTION esat_fn(Ta)
RESULT(esat)
1402 REAL(KIND(1D0))::esat
1404 REAL(KIND(1D0)),
PARAMETER::A = 6.106
1405 REAL(KIND(1D0)),
PARAMETER::B = 17.27
1406 REAL(KIND(1D0)),
PARAMETER::C = 237.3
1408 esat = a*exp(b*ta/(c + ta))
1416 FUNCTION qsat_fn(Ta, pres)
RESULT(qsat)
1419 REAL(KIND(1D0))::qsat
1420 REAL(KIND(1D0))::pres
1422 REAL(KIND(1D0)),
PARAMETER::molar = 0.028965
1423 REAL(KIND(1D0)),
PARAMETER::molar_wat_vap = 0.0180153
1426 qsat = (molar_wat_vap/molar)*es/pres
1433 FUNCTION qa_fn(Ta, RH, pres)
RESULT(qa)
1439 REAL(KIND(1D0))::pres
1442 REAL(KIND(1D0)),
PARAMETER::molar = 0.028965
1443 REAL(KIND(1D0)),
PARAMETER::molar_wat_vap = 0.0180153
1447 qa = (molar_wat_vap/molar)*ea/pres
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...