se_single_column_mod.F90 coverage: 0.00 %func 0.00 %block
1) module se_single_column_mod
2) !--------------------------------------------------------
3) !
4) ! Module for the SE single column model
5)
6) use element_mod, only: element_t
7) use scamMod
8) use constituents, only: cnst_get_ind
9) use dimensions_mod, only: nelemd, np
10) use time_manager, only: get_nstep, dtime, is_first_step
11) use ppgrid, only: begchunk
12) #ifdef MODEL_THETA_L
13) use element_ops, only: get_R_star
14) use eos, only: pnh_and_exner_from_eos
15) #endif
16) use element_ops, only: get_temperature
17)
18) implicit none
19)
20) public scm_setinitial
21) public scm_setfield
22) public apply_SC_forcing
23)
24) !=========================================================================
25) contains
26) !=========================================================================
27)
28) subroutine scm_setinitial(elem)
29)
30) implicit none
31)
32) type(element_t), intent(inout) :: elem(:)
33)
34) integer i, j, k, ie, thelev
35) integer inumliq, inumice, icldliq, icldice
36)
37) if (.not. use_replay .and. get_nstep() .eq. 0) then
38) call cnst_get_ind('NUMLIQ', inumliq, abrtf=.false.)
39) call cnst_get_ind('NUMICE', inumice, abrtf=.false.)
40) call cnst_get_ind('CLDLIQ', icldliq)
41) call cnst_get_ind('CLDICE', icldice)
42)
43) do ie=1,nelemd
44) do j=1,np
45) do i=1,np
46)
47) ! Find level where tobs is no longer zero
48) thelev=1
49) do k=1, PLEV
50) if (tobs(k) .ne. 0) then
51) thelev=k
52) go to 1000
53) endif
54) enddo
55)
56) 1000 continue
57)
58) if (get_nstep() .le. 1) then
59) do k=1,thelev-1
60) #ifdef MODEL_THETA_L
61) tobs(k)=elem(ie)%derived%FT(i,j,k)
62) #else
63) tobs(k)=elem(ie)%state%T(i,j,k,1)
64) #endif
65) qobs(k)=elem(ie)%state%Q(i,j,k,1)
66) enddo
67) else
68) #ifdef MODEL_THETA_L
69) tobs(:)=elem(ie)%derived%FT(i,j,:)
70) #else
71) tobs(:)=elem(ie)%state%T(i,j,:,1)
72) #endif
73) qobs(:)=elem(ie)%state%Q(i,j,:,1)
74) endif
75)
76) if (get_nstep() .eq. 0) then
77) do k=thelev, PLEV
78) #ifdef MODEL_THETA_L
79) if (have_t) elem(ie)%derived%FT(i,j,k)=tobs(k)
80) #else
81) if (have_t) elem(ie)%state%T(i,j,k,1)=tobs(k)
82) #endif
83) if (have_q) elem(ie)%state%Q(i,j,k,1)=qobs(k)
84) enddo
85)
86) do k=1,PLEV
87) if (have_ps) elem(ie)%state%ps_v(i,j,1) = psobs
88) if (have_u) elem(ie)%state%v(i,j,1,k,1) = uobs(k)
89) if (have_v) elem(ie)%state%v(i,j,2,k,1) = vobs(k)
90) if (have_numliq) elem(ie)%state%Q(i,j,k,inumliq) = numliqobs(k)
91) if (have_cldliq) elem(ie)%state%Q(i,j,k,icldliq) = cldliqobs(k)
92) if (have_numice) elem(ie)%state%Q(i,j,k,inumice) = numiceobs(k)
93) if (have_cldice) elem(ie)%state%Q(i,j,k,icldice) = cldiceobs(k)
94) if (have_omega) elem(ie)%derived%omega_p(i,j,k) = wfld(k)
95) enddo
96)
97) endif
98)
99) enddo
100) enddo
101) enddo
102) endif
103)
104) end subroutine scm_setinitial
105)
106) subroutine scm_setfield(elem,iop_update_phase1)
107)
108) implicit none
109)
110) logical, intent(in) :: iop_update_phase1
111) type(element_t), intent(inout) :: elem(:)
112)
113) integer i, j, k, ie
114)
115) do ie=1,nelemd
116) if (have_ps .and. use_replay .and. .not. iop_update_phase1) elem(ie)%state%ps_v(:,:,:) = psobs
117) if (have_ps .and. .not. use_replay) elem(ie)%state%ps_v(:,:,:) = psobs
118) do i=1, PLEV
119) if (have_omega .and. iop_update_phase1) elem(ie)%derived%omega_p(:,:,i)=wfld(i) ! set t to tobs at first
120) end do
121) end do
122)
123) end subroutine scm_setfield
124)
125) subroutine apply_SC_forcing(elem,hvcoord,tl,n,t_before_advance,nets,nete)
126) !
127) use scamMod, only: single_column, use_3dfrc
128) use kinds, only : real_kind
129) use dimensions_mod, only : np, np, nlev, npsq
130) use control_mod, only : use_cpstar, qsplit
131) use hybvcoord_mod, only : hvcoord_t
132) use element_mod, only : element_t
133) use physical_constants, only : Cp, Rgas, cpwater_vapor
134) use time_mod
135) use constituents, only: pcnst
136) use time_manager, only: get_nstep
137) use shr_const_mod, only: SHR_CONST_PI
138)
139) integer :: t1,t2,n,nets,nete,pp
140) type (element_t) , intent(inout), target :: elem(:)
141) type (hvcoord_t) :: hvcoord
142) type (TimeLevel_t), intent(in) :: tl
143) logical :: t_before_advance, do_column_scm
144) real(kind=real_kind), parameter :: rad2deg = 180.0_real_kind / SHR_CONST_PI
145)
146) integer :: ie,k,i,j,t,nm_f
147) real (kind=real_kind), dimension(np,np,nlev) :: dpt1,dpt2 ! delta pressure
148) real (kind=real_kind), dimension(np,np) :: E
149) real (kind=real_kind), dimension(np,np) :: suml,suml2,v1,v2
150) real (kind=real_kind), dimension(np,np,nlev) :: sumlk, suml2k
151) real (kind=real_kind), dimension(np,np,nlev) :: p,T_v,phi, pnh
152) real (kind=real_kind), dimension(np,np,nlev+1) :: dpnh_dp_i
153) real (kind=real_kind), dimension(np,np,nlev) :: dp, exner, vtheta_dp, Rstar
154) real (kind=real_kind), dimension(np,np,nlev) :: dpscm
155) real (kind=real_kind) :: cp_star1,cp_star2,qval_t1,qval_t2
156) real (kind=real_kind) :: Qt,dt
157) real (kind=real_kind), dimension(nlev,pcnst) :: stateQin1, stateQin2, stateQin_qfcst
158) real (kind=real_kind), dimension(nlev,pcnst) :: forecast_q
159) real (kind=real_kind), dimension(nlev) :: temp_tend, dummy2, forecast_t, forecast_u, forecast_v
160) real (kind=real_kind) :: forecast_ps
161) real (kind=real_kind) :: temperature(np,np,nlev)
162) logical :: wet
163)
164) integer:: icount
165)
166) nm_f = 1
167) if (t_before_advance) then
168) t1=tl%nm1
169) t2=tl%n0
170) else
171) t1=tl%n0
172) t2=tl%np1
173) endif
174)
175) ! For SCM only one column is considered
176) ie = 1
177) #if (defined COLUMN_OPENMP)
178) !$omp parallel do private(k)
179) #endif
180)
181) do k=1,nlev
182) p(:,:,k) = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*elem(ie)%state%ps_v(:,:,t1)
183) end do
184)
185) #ifdef MODEL_THETA_L
186) ! If using the theta-l dycore then need to get the exner function
187) ! and reference levels "dp", so we can convert the SCM forecasted
188) ! temperature back potential temperature on reference levels.
189) dp(:,:,:) = elem(ie)%state%dp3d(:,:,:,t1)
190) call pnh_and_exner_from_eos(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,t1),&
191) dp,elem(ie)%state%phinh_i(:,:,:,t1),pnh,exner,dpnh_dp_i)
192) #endif
193)
194) dt=dtime
195)
196) ! Get temperature from dynamics state
197) call get_temperature(elem(ie),temperature,hvcoord,t1)
198)
199) ! For SCM we only consider one point
200) i=1
201) j=1
202)
203) stateQin_qfcst(:,:) = elem(ie)%state%Q(i,j,:,:)
204) stateQin1(:,:) = stateQin_qfcst(:,:)
205) stateQin2(:,:) = stateQin_qfcst(:,:)
206)
207) if (.not. use_3dfrc) then
208) temp_tend(:) = 0.0_real_kind
209) else
210) temp_tend(:) = elem(ie)%derived%fT(i,j,:)
211) endif
212) dummy2(:) = 0.0_real_kind
213) forecast_ps = elem(ie)%state%ps_v(i,j,t1)
214)
215) #ifdef MODEL_THETA_L
216) ! At first time step the tendency term is set to the
217) ! initial condition temperature profile. Do NOT use
218) ! as it will cause temperature profile to blow up.
219) if ( is_first_step() ) then
220) temp_tend(:) = 0.0
221) endif
222) #endif
223)
224) call forecast(begchunk,elem(ie)%state%ps_v(i,j,t1),&
225) elem(ie)%state%ps_v(i,j,t1),forecast_ps,forecast_u,&
226) elem(ie)%state%v(i,j,1,:,t1),elem(ie)%state%v(i,j,1,:,t1),&
227) forecast_v,elem(ie)%state%v(i,j,2,:,t1),&
228) elem(ie)%state%v(i,j,2,:,t1),forecast_t,&
229) temperature(i,j,:),temperature(i,j,:),&
230) forecast_q,stateQin2,stateQin1,dt,temp_tend,dummy2,dummy2,&
231) stateQin_qfcst,p(i,j,:),stateQin1,1)
232)
233) elem(ie)%state%Q(i,j,:,:) = forecast_q(:,:)
234)
235) #ifdef MODEL_THETA_L
236) ! If running theta-l model then the forecast temperature needs
237) ! to be converted back to potential temperature on reference levels,
238) ! which is what dp_coupling expects
239) call get_R_star(Rstar,elem(ie)%state%Q(:,:,:,1))
240) elem(ie)%state%vtheta_dp(i,j,:,t1) = (forecast_t(:)*Rstar(i,j,:)*dp(i,j,:))/&
241) (Rgas*exner(i,j,:))
242) #else
243) elem(ie)%state%T(i,j,:,t1) = forecast_t(:)
244) #endif
245) elem(ie)%state%v(i,j,1,:,t1) = forecast_u(:)
246) elem(ie)%state%v(i,j,2,:,t1) = forecast_v(:)
247)
248) end subroutine apply_SC_forcing
249)
250) end module se_single_column_mod
