Line data Source code
1 : module apply_iop_forcing_mod
2 :
3 : use shr_kind_mod, only:r8 => shr_kind_r8, i8 => shr_kind_i8
4 : use pmgrid, only:plev, plevp, plon
5 : use constituents, only:pcnst, cnst_get_ind, cnst_name
6 : use physconst, only:rair,cpair
7 : use cam_logfile, only:iulog
8 : use hybvcoord_mod, only: hvcoord_t
9 : use scamMod, only: use_3dfrc, single_column, have_u, have_v, divT3d, divq3d, divt, divq, &
10 : wfld, uobs, vobs, tobs, qobs, plevs0, have_divt3d, have_divq3d, &
11 : scm_relax_bot_p,scm_relax_linear,scm_relax_tau_bot_sec, &
12 : scm_relax_tau_sec,scm_relax_tau_top_sec,scm_relax_top_p, &
13 : scm_relaxation,scm_relax_fincl,qinitobs
14 :
15 : use cam_abortutils, only: endrun
16 : use string_utils, only: to_upper
17 :
18 : implicit none
19 :
20 : public advance_iop_forcing
21 : public advance_iop_nudging
22 :
23 : !=========================================================================
24 : contains
25 : !=========================================================================
26 :
27 0 : subroutine advance_iop_forcing(scm_dt, ps_in, & ! In
28 : u_in, v_in, t_in, q_in, t_phys_frc, q_phys_frc, hvcoord, & ! In
29 : u_update, v_update, t_update, q_update) ! Out
30 :
31 : !-----------------------------------------------------------------------
32 : !
33 : ! Purpose:
34 : ! Apply large scale forcing for t, q, u, and v as provided by the
35 : ! case IOP forcing file.
36 : !
37 : ! Author:
38 : ! Original version: Adopted from CAM3.5/CAM5
39 : ! Updated version for E3SM: Peter Bogenschutz (bogenschutz1@llnl.gov)
40 : ! and replaces the forecast.F90 routine in CAM3.5/CAM5/CAM6/E3SMv1/E3SMv2
41 : !
42 : !-----------------------------------------------------------------------
43 :
44 : ! Input arguments
45 : real(r8), intent(in) :: ps_in ! surface pressure [Pa]
46 : real(r8), intent(in) :: u_in(plev) ! zonal wind [m/s]
47 : real(r8), intent(in) :: v_in(plev) ! meridional wind [m/s]
48 : real(r8), intent(in) :: t_in(plev) ! temperature [K]
49 : real(r8), intent(in) :: q_in(plev,pcnst) ! q tracer array [units vary] already vertically advected
50 : real(r8), intent(in) :: t_phys_frc(plev) ! temperature forcing from physics [K/s]
51 : real(r8), intent(in) :: q_phys_frc(plev,pcnst) ! change in q due to physics.
52 : type (hvcoord_t), intent(in) :: hvcoord
53 : real(r8), intent(in) :: scm_dt ! model time step [s]
54 :
55 : ! Output arguments
56 : real(r8), intent(out) :: t_update(plev) ! updated temperature [K]
57 : real(r8), intent(out) :: q_update(plev,pcnst)! updated q tracer array [units vary]
58 : real(r8), intent(out) :: u_update(plev) ! updated zonal wind [m/s]
59 : real(r8), intent(out) :: v_update(plev) ! updated meridional wind [m/s]
60 :
61 : ! Local variables
62 : real(r8) pmidm1(plev) ! pressure at model levels
63 : real(r8) pintm1(plevp) ! pressure at model interfaces
64 : real(r8) pdelm1(plev) ! pdel(k) = pint (k+1)-pint (k)
65 : real(r8) t_lsf(plev) ! storage for temperature large scale forcing
66 : real(r8) q_lsf(plev,pcnst) ! storage for moisture large scale forcing
67 : real(r8) fac, t_expan
68 :
69 : integer i,k,m ! longitude, level, constituent indices
70 :
71 : character(len=*), parameter :: subname = 'advance_iop_forcing'
72 :
73 : ! Get vertical level profiles
74 0 : call plevs0(plev, ps_in, hvcoord%ps0, hvcoord%hyam, hvcoord%hybm, hvcoord%hyai, hvcoord%hybi, pintm1 ,pmidm1 ,pdelm1)
75 :
76 : ! Advance T and Q due to large scale forcing
77 0 : if (use_3dfrc) then
78 0 : if(.not.(have_divt3d.and.have_divq3d)) call endrun(subname//": FATAL: divt3d and divq3d not available")
79 0 : t_lsf(:plev) = divt3d(:plev)
80 0 : q_lsf(:plev,:pcnst) = divq3d(:plev,:pcnst)
81 : else
82 0 : t_lsf(:plev) = divt(:plev)
83 0 : q_lsf(:plev,:pcnst) = divq(:plev,:pcnst)
84 : endif
85 :
86 0 : do k=1,plev
87 : ! Initialize thermal expansion term to zero. This term is only
88 : ! considered if three dimensional forcing is not provided by IOP forcing file.
89 0 : t_expan = 0._r8
90 :
91 0 : if (.not. use_3dfrc) then
92 0 : t_expan = scm_dt*wfld(k)*t_in(k)*rair/(cpair*pmidm1(k))
93 : endif
94 :
95 0 : if (use_3dfrc) then
96 0 : do m=1,pcnst
97 : ! When using 3d dynamics tendencies, SCM skips the vertical advection step and thus
98 : ! q_in at this point has not had physics tendencies applied
99 0 : q_update(k,m) = q_in(k,m) + scm_dt*(q_phys_frc(k,m) + q_lsf(k,m))
100 : end do
101 0 : t_update(k) = t_in(k) + t_expan + scm_dt*(t_phys_frc(k) + t_lsf(k))
102 : else
103 0 : do m=1,pcnst
104 : ! When not using 3d dynamics tendencies, q_in at this point has had physics tend
105 : ! applied and has been vertically advected. Only horizontal dyn tend needed for forecast.
106 0 : q_update(k,m) = q_in(k,m) + scm_dt*q_lsf(k,m)
107 : end do
108 0 : t_update(k) = t_in(k) + t_expan + scm_dt*t_lsf(k)
109 : end if
110 : end do
111 :
112 : ! Set U and V fields
113 :
114 0 : if ( have_v .and. have_u ) then
115 0 : do k=1,plev
116 0 : u_update(k) = uobs(k)
117 0 : v_update(k) = vobs(k)
118 : enddo
119 : endif
120 :
121 0 : end subroutine advance_iop_forcing
122 :
123 : !=========================================================================
124 :
125 0 : subroutine advance_iop_nudging(ztodt, ps_in, & ! In
126 : tfcst, qfcst, ufcst, vfcst, hvcoord, & ! Inout
127 : relaxt, relaxq ) ! Out
128 :
129 : !-----------------------------------------------------------------------
130 : !
131 : ! Purpose:
132 : ! Option to nudge t and q to observations as specified by the IOP file
133 : !-----------------------------------------------------------------------
134 :
135 : ! Input arguments
136 : real(r8), intent(in) :: ztodt ! model time step [s]
137 : real(r8), intent(in) :: ps_in ! surface pressure [Pa]
138 : type (hvcoord_t), intent(in) :: hvcoord
139 :
140 : ! Output arguments
141 : real(r8), intent(inout) :: tfcst(plev) ! updated temperature [K]
142 : real(r8), intent(inout) :: qfcst(plon,plev,pcnst) ! updated const field
143 : real(r8), intent(inout) :: ufcst(plev) ! updated U wind
144 : real(r8), intent(inout) :: vfcst(plev) ! updated V wind
145 : real(r8), intent(out) :: relaxt(plev) ! relaxation of temperature [K/s]
146 : real(r8), intent(out) :: relaxq(plev) ! relaxation of vapor [kg/kg/s]
147 :
148 : ! Local variables
149 : integer :: i, k, m
150 : real(r8) pmidm1(plev) ! pressure at model levels
151 : real(r8) pintm1(plevp) ! pressure at model interfaces
152 : real(r8) pdelm1(plev) ! pdel(k) = pint (k+1)-pint (k)
153 :
154 : ! --------------------------- !
155 : ! For 'scm_relaxation' switch !
156 : ! --------------------------- !
157 :
158 : real(r8) rtau(plev)
159 : real(r8) relax_T(plev)
160 : real(r8) relax_u(plev)
161 : real(r8) relax_v(plev)
162 : real(r8) relax_q(plev,pcnst)
163 : ! +++BPM: allow linear relaxation profile
164 : real(r8) rslope ! [optional] slope for linear relaxation profile
165 : real(r8) rycept ! [optional] y-intercept for linear relaxtion profile
166 : logical scm_fincl_empty
167 :
168 : ! ------------------------------------------------------------------- !
169 : ! Relaxation to the observed or specified state !
170 : ! We should specify relaxation time scale ( rtau ) and !
171 : ! target-relaxation state ( in the current case, either 'obs' or 0 ) !
172 : ! ------------------------------------------------------------------- !
173 :
174 0 : if ( .not. scm_relaxation) return
175 :
176 0 : call plevs0(plev, ps_in, hvcoord%ps0, hvcoord%hyam, hvcoord%hybm, hvcoord%hyai, hvcoord%hybi, pintm1 ,pmidm1 ,pdelm1)
177 :
178 0 : relax_T(:) = 0._r8
179 0 : relax_u(:) = 0._r8
180 0 : relax_v(:) = 0._r8
181 0 : relax_q(:plev,:pcnst) = 0._r8
182 : ! +++BPM: allow linear relaxation profile
183 : ! scm_relaxation is a logical from scamMod
184 : ! scm_relax_tau_top_sec and scm_relax_tau_bot_sec are the relaxation times at top and bottom of layer
185 : ! also defined in scamMod
186 0 : if ( scm_relax_linear ) then
187 0 : rslope = (scm_relax_top_p - scm_relax_bot_p)/(scm_relax_tau_top_sec - scm_relax_tau_bot_sec)
188 0 : rycept = scm_relax_tau_top_sec - (rslope*scm_relax_top_p)
189 : endif
190 :
191 0 : scm_fincl_empty=.true.
192 0 : do i=1,pcnst
193 0 : if (len_trim(scm_relax_fincl(i)) > 0) then
194 0 : scm_fincl_empty=.false.
195 0 : scm_relax_fincl(i)=trim(to_upper(scm_relax_fincl(i)))
196 : end if
197 : end do
198 :
199 0 : do k = 1, plev
200 0 : if ( pmidm1(k) <= scm_relax_bot_p.and.pmidm1(k) >= scm_relax_top_p ) then ! inside layer
201 0 : if (scm_relax_linear) then
202 0 : rtau(k) = rslope*pmidm1(k) + rycept ! linear regime
203 : else
204 0 : rtau(k) = max( ztodt, scm_relax_tau_sec ) ! constant for whole layer / no relax outside
205 : endif
206 0 : else if (scm_relax_linear .and. pmidm1(k) <= scm_relax_top_p ) then ! not linear => do nothing / linear => use upper value
207 0 : rtau(k) = scm_relax_tau_top_sec ! above layer keep rtau equal to the top
208 : endif
209 : ! +BPM: this can't be the best way...
210 : ! I put this in because if rtau doesn't get set above, then I don't want to do any relaxation in that layer.
211 : ! maybe the logic of this whole loop needs to be re-thinked.
212 0 : if (rtau(k) /= 0) then
213 0 : relax_T(k) = - ( tfcst(k) - tobs(k) ) / rtau(k)
214 0 : relax_u(k) = - ( ufcst(k) - uobs(k) ) / rtau(k)
215 0 : relax_v(k) = - ( vfcst(k) - vobs(k) ) / rtau(k)
216 0 : relax_q(k,1) = - ( qfcst(1,k,1) - qobs(k) ) / rtau(k)
217 0 : do m = 2, pcnst
218 0 : relax_q(k,m) = - ( qfcst(1,k,m) - qinitobs(k,m) ) / rtau(k)
219 : enddo
220 0 : if (scm_fincl_empty .or. ANY(scm_relax_fincl(:) == 'T')) &
221 0 : tfcst(k) = tfcst(k) + relax_T(k) * ztodt
222 0 : if (scm_fincl_empty .or.ANY(scm_relax_fincl(:) == 'U')) &
223 0 : ufcst(k) = ufcst(k) + relax_u(k) * ztodt
224 0 : if (scm_fincl_empty .or. ANY(scm_relax_fincl(:) == 'V')) &
225 0 : vfcst(k) = vfcst(k) + relax_v(k) * ztodt
226 0 : do m = 1, pcnst
227 0 : if (scm_fincl_empty .or. ANY(scm_relax_fincl(:) == trim(to_upper(cnst_name(m)))) ) then
228 0 : qfcst(1,k,m) = qfcst(1,k,m) + relax_q(k,m) * ztodt
229 : end if
230 : enddo
231 : end if
232 : enddo
233 :
234 : end subroutine advance_iop_nudging
235 :
236 : !-----------------------------------------------------------------------
237 :
238 : end module apply_iop_forcing_mod
|
