Line data Source code
1 : module se_single_column_mod
2 : !--------------------------------------------------------
3 : !
4 : ! Module for the SE single column model
5 :
6 : use shr_kind_mod, only: r8=>shr_kind_r8
7 : use element_mod, only: element_t
8 : use scamMod, only: have_t, have_q, have_u, have_v, have_ps, have_numliq, &
9 : have_cldliq, have_numice, have_cldice, have_omega, use_camiop, &
10 : tobs, qobs,have_numliq, numliqobs, cldliqobs, numiceobs, cldiceobs, &
11 : wfld, psobs,uobs,vobs,tobs,divt,divQ,divT3d,divq3d,precobs,lhflxobs, &
12 : shflxobs, tground, have_ps, have_tg, have_lhflx, have_shflx, have_t, &
13 : have_omega, have_cldliq, have_divt, have_divq, have_divt3d, have_divq3d, &
14 : use_3dfrc,scmlat,scmlon
15 : use constituents, only: cnst_get_ind, pcnst
16 : use dimensions_mod, only: nelemd, np, nlev, qsize
17 : use time_manager, only: get_nstep, is_first_step, get_step_size, is_first_restart_step
18 : use ppgrid, only: begchunk
19 : use se_dyn_time_mod, only: timelevel_qdp
20 : use cam_history, only: outfld
21 :
22 : implicit none
23 :
24 : private
25 : save
26 :
27 : public scm_setinitial
28 : public scm_setfield
29 : public apply_SC_forcing
30 : public iop_broadcast
31 : public scm_dyn_grid_indicies
32 :
33 : integer, public :: indx_scm, ie_scm, i_scm, j_scm
34 :
35 : integer :: tl_f, tl_fqdp, thelev
36 :
37 : !=========================================================================
38 : contains
39 : !=========================================================================
40 :
41 0 : subroutine scm_setinitial(elem)
42 :
43 : use dyn_grid, only: TimeLevel
44 : use control_mod, only: qsplit
45 :
46 : implicit none
47 :
48 : type(element_t), intent(inout) :: elem(:)
49 :
50 : integer :: k
51 : integer :: inumliq, inumice, icldliq, icldice
52 :
53 0 : call scm_dyn_grid_indicies(elem,scmlat,scmlon,ie_scm,i_scm,j_scm,indx_scm)
54 :
55 0 : tl_f = timelevel%n0
56 0 : call TimeLevel_Qdp(timelevel, qsplit, tl_fqdp)
57 :
58 0 : if (.not. use_camiop .and. get_nstep() == 0) then
59 0 : call cnst_get_ind('NUMLIQ', inumliq, abort=.false.)
60 0 : call cnst_get_ind('NUMICE', inumice, abort=.false.)
61 0 : call cnst_get_ind('CLDLIQ', icldliq)
62 0 : call cnst_get_ind('CLDICE', icldice)
63 :
64 : ! Find level where tobs is no longer zero
65 0 : thelev=minloc(abs(tobs), 1, mask=abs(tobs) > 0)
66 :
67 0 : if (get_nstep() <= 1) then
68 0 : do k=1,thelev-1
69 0 : tobs(k)=elem(ie_scm)%state%T(i_scm,j_scm,k,tl_f)
70 0 : qobs(k)=elem(ie_scm)%state%qdp(i_scm,j_scm,k,1,tl_fqdp)/elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
71 : enddo
72 : else
73 0 : tobs(:)=elem(ie_scm)%state%T(i_scm,j_scm,:,tl_f)
74 0 : qobs(:)=elem(ie_scm)%state%qdp(i_scm,j_scm,:,1,tl_fqdp)/elem(ie_scm)%state%dp3d(i_scm,j_scm,:,tl_f)
75 : endif
76 :
77 0 : if (get_nstep() == 0) then
78 0 : do k=thelev, NLEV
79 0 : if (have_t) elem(ie_scm)%state%T(i_scm,j_scm,k,tl_f)=tobs(k)
80 0 : if (have_q) elem(ie_scm)%state%qdp(i_scm,j_scm,k,1,tl_fqdp)=qobs(k)*elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
81 : enddo
82 :
83 0 : do k=1,NLEV
84 0 : if (have_ps) elem(ie_scm)%state%psdry(i_scm,j_scm) = psobs
85 0 : if (have_u) elem(ie_scm)%state%v(i_scm,j_scm,1,k,tl_f) = uobs(k)
86 0 : if (have_v) elem(ie_scm)%state%v(i_scm,j_scm,2,k,tl_f) = vobs(k)
87 0 : if (have_numliq) elem(ie_scm)%state%qdp(i_scm,j_scm,k,inumliq,tl_fqdp) = &
88 0 : numliqobs(k)*elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
89 0 : if (have_cldliq) elem(ie_scm)%state%qdp(i_scm,j_scm,k,icldliq,tl_fqdp) = &
90 0 : cldliqobs(k)*elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
91 0 : if (have_numice) elem(ie_scm)%state%qdp(i_scm,j_scm,k,inumice,tl_fqdp) = &
92 0 : numiceobs(k)*elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
93 0 : if (have_cldice) elem(ie_scm)%state%qdp(i_scm,j_scm,k,icldice,tl_fqdp) = &
94 0 : cldiceobs(k)*elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
95 0 : if (have_omega) elem(ie_scm)%derived%omega(i_scm,j_scm,k) = wfld(k)
96 : enddo
97 :
98 : endif
99 :
100 : endif
101 :
102 0 : end subroutine scm_setinitial
103 :
104 0 : subroutine scm_setfield(elem,iop_update_phase1)
105 :
106 : !---------------------------------------------------------
107 : ! Purpose: Update various fields based on available data
108 : ! provided by IOP file
109 : !----------------------------------------------------------
110 :
111 0 : use control_mod, only: qsplit
112 : use dyn_grid, only: TimeLevel
113 :
114 : implicit none
115 :
116 : logical, intent(in) :: iop_update_phase1
117 : type(element_t), intent(inout) :: elem(:)
118 :
119 : integer :: k
120 : integer :: tl_f, tl_fqdp
121 :
122 0 : tl_f = timelevel%n0
123 0 : call TimeLevel_Qdp(timelevel, qsplit, tl_fqdp)
124 :
125 0 : if (have_ps .and. use_camiop .and. .not. iop_update_phase1) elem(ie_scm)%state%psdry(:,:) = psobs
126 0 : if (have_ps .and. .not. use_camiop) elem(ie_scm)%state%psdry(:,:) = psobs
127 0 : do k=1, NLEV
128 0 : if (have_omega .and. iop_update_phase1) elem(ie_scm)%derived%omega(:,:,k)=wfld(k) ! set t to tobs at first
129 0 : if (k < thelev) then
130 0 : tobs(k) = elem(ie_scm)%state%T(i_scm,j_scm,k,tl_f)
131 0 : qobs(k) = elem(ie_scm)%state%qdp(i_scm,j_scm,k,1,tl_fqdp)/elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f)
132 0 : uobs(k) = elem(ie_scm)%state%v(i_scm,j_scm,1,k,tl_f)
133 0 : vobs(k) = elem(ie_scm)%state%v(i_scm,j_scm,2,k,tl_f)
134 : end if
135 : end do
136 :
137 0 : end subroutine scm_setfield
138 :
139 0 : subroutine apply_SC_forcing(elem,hvcoord,tl,n,t_before_advance)
140 : !
141 0 : use scamMod, only: single_column, use_3dfrc
142 : use hybvcoord_mod, only: hvcoord_t
143 : use se_dyn_time_mod,only: TimeLevel_t
144 : use control_mod, only: qsplit
145 : use apply_iop_forcing_mod, only:advance_iop_forcing, advance_iop_nudging
146 :
147 : type (element_t), intent(inout), target :: elem(:)
148 : type (hvcoord_t), intent(in) :: hvcoord
149 : type (TimeLevel_t), intent(in) :: tl
150 : logical, intent(in) :: t_before_advance
151 : integer, intent(in) :: n
152 :
153 : integer :: k, m
154 : real (r8) :: dt
155 : logical :: iop_nudge_tq = .false.
156 : real (r8), dimension(nlev,pcnst) :: stateQ_in, q_update, q_phys_frc
157 : real (r8), dimension(nlev) :: t_phys_frc, t_update, u_update, v_update
158 : real (r8), dimension(nlev) :: t_in, u_in, v_in
159 : real (r8), dimension(nlev) :: relaxt, relaxq
160 : real (r8), dimension(nlev) :: tdiff_dyn, qdiff_dyn
161 :
162 : !-----------------------------------------------------------------------
163 :
164 0 : tl_f = tl%n0
165 :
166 0 : call TimeLevel_Qdp(tl, qsplit, tl_fqdp)
167 :
168 0 : dt = get_step_size()
169 :
170 : ! Set initial profiles for current column
171 0 : do m=1,pcnst
172 0 : stateQ_in(:nlev,m) = elem(ie_scm)%state%Qdp(i_scm,j_scm,:nlev,m,tl_fqdp)/elem(ie_scm)%state%dp3d(i_scm,j_scm,:nlev,tl_f)
173 : end do
174 0 : t_in(:nlev) = elem(ie_scm)%state%T(i_scm,j_scm,:nlev,tl_f)
175 0 : u_in(:nlev) = elem(ie_scm)%state%v(i_scm,j_scm,1,:nlev,tl_f)
176 0 : v_in(:nlev) = elem(ie_scm)%state%v(i_scm,j_scm,2,:nlev,tl_f)
177 :
178 0 : t_phys_frc(:) = elem(ie_scm)%derived%fT(i_scm,j_scm,:)
179 0 : q_phys_frc(:,:qsize) = elem(ie_scm)%derived%fQ(i_scm,j_scm,:,:qsize)/dt
180 :
181 : ! Call the main subroutine to update t, q, u, and v according to
182 : ! large scale forcing as specified in IOP file.
183 : call advance_iop_forcing(dt,elem(ie_scm)%state%psdry(i_scm,j_scm),& ! In
184 : u_in,v_in,t_in,stateQ_in,t_phys_frc, q_phys_frc, hvcoord, & ! In
185 0 : u_update,v_update,t_update,q_update) ! Out
186 :
187 : ! Nudge to observations if desired, for T & Q only if in SCM mode
188 0 : if (iop_nudge_tq ) then
189 0 : call advance_iop_nudging(dt,elem(ie_scm)%state%psdry(i_scm,j_scm),& ! In
190 : t_update,q_update,u_update,v_update, hvcoord, & ! Inout
191 0 : relaxt,relaxq) ! Out
192 : endif
193 :
194 0 : if (use_3dfrc) then ! vertical remap of dynamics not run need to update state%dp3d using new psdry
195 0 : do k=1,nlev
196 0 : elem(ie_scm)%state%dp3d(i_scm,j_scm,k,tl_f) = (hvcoord%hyai(k+1)-hvcoord%hyai(k))*hvcoord%ps0 + &
197 0 : (hvcoord%hybi(k+1)-hvcoord%hybi(k))*elem(ie_scm)%state%psdry(i_scm,j_scm)
198 : end do
199 : end if
200 :
201 : ! Update qdp using new dp3d
202 0 : do m=1,pcnst
203 : ! Update the Qdp array
204 0 : elem(ie_scm)%state%Qdp(i_scm,j_scm,:nlev,m,tl_fqdp) = &
205 0 : q_update(:nlev,m) * elem(ie_scm)%state%dp3d(i_scm,j_scm,:nlev,tl_f)
206 : enddo
207 :
208 : ! Update prognostic variables to the current values
209 0 : elem(ie_scm)%state%T(i_scm,j_scm,:,tl_f) = t_update(:)
210 0 : elem(ie_scm)%state%v(i_scm,j_scm,1,:,tl_f) = u_update(:)
211 0 : elem(ie_scm)%state%v(i_scm,j_scm,2,:,tl_f) = v_update(:)
212 :
213 : ! Evaluate the differences in state information from observed
214 : ! (done for diganostic purposes only)
215 0 : do k = 1, nlev
216 0 : tdiff_dyn(k) = t_update(k) - tobs(k)
217 0 : qdiff_dyn(k) = q_update(k,1) - qobs(k)
218 : end do
219 :
220 : ! Add various diganostic outfld calls
221 0 : call outfld('TDIFF',tdiff_dyn,1,begchunk)
222 0 : call outfld('QDIFF',qdiff_dyn,1,begchunk)
223 0 : call outfld('TOBS',tobs,1,begchunk)
224 0 : call outfld('QOBS',qobs,1,begchunk)
225 0 : call outfld('DIVQ',divq,1,begchunk)
226 0 : call outfld('DIVT',divt,1,begchunk)
227 0 : call outfld('DIVQ3D',divq3d,1,begchunk)
228 0 : call outfld('DIVT3D',divt3d,1,begchunk)
229 0 : call outfld('PRECOBS',precobs,1,begchunk)
230 0 : call outfld('LHFLXOBS',lhflxobs,1,begchunk)
231 0 : call outfld('SHFLXOBS',shflxobs,1,begchunk)
232 :
233 0 : call outfld('TRELAX',relaxt,1,begchunk)
234 0 : call outfld('QRELAX',relaxq,1,begchunk)
235 :
236 :
237 0 : end subroutine apply_SC_forcing
238 : !=========================================================================
239 0 : subroutine iop_broadcast()
240 :
241 : !---------------------------------------------------------
242 : ! Purpose: Broadcast relevant logical
243 : ! flags and data to all processors
244 : !----------------------------------------------------------
245 :
246 0 : use spmd_utils, only: mpi_logical, mpi_real8, masterproc, iam, mpicom, mstrid=>masterprocid
247 : use cam_abortutils, only: endrun
248 :
249 : integer :: ierr
250 : character(len=*), parameter :: sub = 'radiation_readnl'
251 :
252 : #ifdef SPMD
253 0 : call mpi_bcast(have_ps,1,mpi_logical,mstrid,mpicom,ierr)
254 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_ps")
255 0 : call mpi_bcast(have_tg,1,mpi_logical,mstrid,mpicom,ierr)
256 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_tg")
257 0 : call mpi_bcast(have_lhflx,1,mpi_logical,mstrid,mpicom,ierr)
258 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_lhflx")
259 0 : call mpi_bcast(have_shflx,1,mpi_logical,mstrid,mpicom,ierr)
260 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_shflx")
261 0 : call mpi_bcast(have_t,1,mpi_logical,mstrid,mpicom,ierr)
262 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_t")
263 0 : call mpi_bcast(have_q,1,mpi_logical,mstrid,mpicom,ierr)
264 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_q")
265 0 : call mpi_bcast(have_u,1,mpi_logical,mstrid,mpicom,ierr)
266 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_u")
267 0 : call mpi_bcast(have_v,1,mpi_logical,mstrid,mpicom,ierr)
268 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_v")
269 0 : call mpi_bcast(have_omega,1,mpi_logical,mstrid,mpicom,ierr)
270 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_omega")
271 0 : call mpi_bcast(have_cldliq,1,mpi_logical,mstrid,mpicom,ierr)
272 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_cldliq")
273 0 : call mpi_bcast(have_divt,1,mpi_logical,mstrid,mpicom,ierr)
274 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_divt")
275 0 : call mpi_bcast(have_divq,1,mpi_logical,mstrid,mpicom,ierr)
276 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_divq")
277 0 : call mpi_bcast(have_divt3d,1,mpi_logical,mstrid,mpicom,ierr)
278 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_divt3d")
279 0 : call mpi_bcast(have_divq3d,1,mpi_logical,mstrid,mpicom,ierr)
280 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: have_divq3d")
281 0 : call mpi_bcast(use_3dfrc,1,mpi_logical,mstrid,mpicom,ierr)
282 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_3dfrc")
283 :
284 0 : call mpi_bcast(psobs,1,mpi_real8,mstrid,mpicom,ierr)
285 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: psobs")
286 0 : call mpi_bcast(tground,1,mpi_real8,mstrid,mpicom,ierr)
287 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: tground")
288 0 : call mpi_bcast(lhflxobs,1,mpi_real8,mstrid,mpicom,ierr)
289 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: lhflxobs")
290 0 : call mpi_bcast(shflxobs,1,mpi_real8,mstrid,mpicom,ierr)
291 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: shflxobs")
292 :
293 0 : call mpi_bcast(tobs,nlev,mpi_real8,mstrid,mpicom,ierr)
294 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: tobs")
295 0 : call mpi_bcast(qobs,nlev,mpi_real8,mstrid,mpicom,ierr)
296 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: qobs")
297 0 : call mpi_bcast(uobs,nlev,mpi_real8,mstrid,mpicom,ierr)
298 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: uobs")
299 0 : call mpi_bcast(vobs,nlev,mpi_real8,mstrid,mpicom,ierr)
300 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: vobs")
301 0 : call mpi_bcast(cldliqobs,nlev,mpi_real8,mstrid,mpicom,ierr)
302 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: cldliqobs")
303 0 : call mpi_bcast(wfld,nlev,mpi_real8,mstrid,mpicom,ierr)
304 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: wfld")
305 :
306 0 : call mpi_bcast(divt,nlev,mpi_real8,mstrid,mpicom,ierr)
307 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: divt")
308 0 : call mpi_bcast(divq,nlev,mpi_real8,mstrid,mpicom,ierr)
309 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: divq")
310 0 : call mpi_bcast(divt3d,nlev,mpi_real8,mstrid,mpicom,ierr)
311 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: divt3d")
312 0 : call mpi_bcast(divq3d,nlev,mpi_real8,mstrid,mpicom,ierr)
313 0 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: divq3d")
314 :
315 : #endif
316 :
317 0 : end subroutine iop_broadcast
318 :
319 : !=========================================================================
320 0 : subroutine scm_dyn_grid_indicies(elem,scmlat,scmlon,ie_scm,i_scm,j_scm,indx_scm)
321 :
322 : !---------------------------------------------------------
323 : ! Purpose: Determine closest column index in the IOP file
324 : ! based on the input scm latitude and longitude
325 : !----------------------------------------------------------
326 :
327 : use shr_const_mod, only: SHR_CONST_PI
328 : use cam_abortutils, only: endrun
329 :
330 : type(element_t), intent(in) :: elem(:)
331 : real (r8), intent(in) :: scmlat,scmlon
332 : integer, intent(out) :: ie_scm, j_scm, i_scm, indx_scm
333 :
334 : integer :: i, j, indx, ie
335 : real(r8) :: scmposlon, minpoint, testlat, testlon, testval
336 : integer :: ierr
337 : real(r8), parameter :: rad2deg = 180.0_r8 / SHR_CONST_PI
338 : character(len=*), parameter :: sub = 'scm_dyn_grid_indicies'
339 :
340 0 : ie_scm=0
341 0 : i_scm=0
342 0 : j_scm=0
343 0 : indx_scm=0
344 0 : minpoint = 1000
345 0 : scmposlon = mod(scmlon + 360._r8,360._r8)
346 0 : do ie=1, nelemd
347 : indx=1
348 0 : do j=1, np
349 0 : do i=1, np
350 0 : testlat=elem(ie)%spherep(i,j)%lat * rad2deg
351 0 : testlon=elem(ie)%spherep(i,j)%lon * rad2deg
352 0 : if (testlon < 0._r8) testlon=testlon+360._r8
353 0 : testval=abs(scmlat-testlat)+abs(scmposlon-testlon)
354 0 : if (testval < minpoint) then
355 0 : ie_scm=ie
356 0 : indx_scm=indx
357 0 : i_scm=i
358 0 : j_scm=j
359 0 : minpoint=testval
360 0 : if (minpoint < 1.e-7_r8) minpoint=0._r8
361 : endif
362 0 : indx=indx+1
363 : enddo
364 : enddo
365 : enddo
366 :
367 0 : if (ie_scm == 0 .or. i_scm == 0 .or. j_scm == 0 .or. indx_scm == 0) then
368 0 : call endrun(sub//':FATAL: Could not find closest SCM point on input datafile')
369 : endif
370 :
371 0 : end subroutine scm_dyn_grid_indicies
372 :
373 : end module se_single_column_mod
|
