FV3 Bundle
fv3jedi_lm_turbulence_mod.F90
Go to the documentation of this file.
2 
6 
7 ! BL Schemes
8 use bldriver
9 
10 ! Saturation table
11 use qsat_util
12 
15 
16 !> Turbulence module
17 
18 implicit none
19 private
21 
22 !> Local trajectory object
24  real(kind_real), allocatable, dimension(:,:,:) :: akq, bkq, ckq
25  real(kind_real), allocatable, dimension(:,:,:) :: aks, bks, cks
26  real(kind_real), allocatable, dimension(:,:,:) :: akv, bkv, ckv
27  real(kind_real), allocatable, dimension(:,:,:) :: pk
28  logical :: set = .false.
30 
31 !> Local constants object
33  real(kind_real), dimension(22) :: turbparams
34  integer, dimension(4) :: turbparamsi
36 
37 !> Turbulence class (self)
39  type(local_traj_turbulence), allocatable :: ltraj(:)
40  type(local_cnst_turbulence) :: lcnst
41  contains
42  procedure :: create
43  procedure :: init_nl
44  procedure :: init_tl
45  procedure :: init_ad
46  procedure :: step_nl
47  procedure :: step_tl
48  procedure :: step_ad
49  procedure :: delete
51 
52 contains
53 
54 ! ------------------------------------------------------------------------------
55 
56 subroutine create(self,conf)
57 
58  implicit none
59 
60  class(fv3jedi_lm_turbulence_type), target, intent(inout) :: self
61  type(fv3jedi_lm_conf), intent(in) :: conf
62 
63  real(kind_real), allocatable, dimension(:) :: pref
64  integer :: l, n
65 
66  if (conf%saveltraj) then
67  allocate(self%ltraj(conf%nt))
68  do n = 1,conf%nt
69  call allocate_ltraj(conf%im,conf%jm,conf%lm,self%ltraj(n))
70  enddo
71  else
72  allocate(self%ltraj(1))
73  call allocate_ltraj(conf%im,conf%jm,conf%lm,self%ltraj(1))
74  endif
75 
76  allocate(pref(0:conf%lm))
77  DO l = 0,conf%lm
78  pref(l) = conf%AK(l+1) + conf%BK(l+1)*p00
79  enddo
80 
81  !Turbulence Parameters
82  self%lcnst%TurbParams(1) = 5.0_kind_real !LOUIS
83  self%lcnst%TurbParams(2) = 160.0_kind_real !LAMBDAM
84  self%lcnst%TurbParams(3) = 1.0_kind_real !LAMBDAM2
85  self%lcnst%TurbParams(4) = 160.0_kind_real !LAMBDAH
86  self%lcnst%TurbParams(5) = 1.0_kind_real !LAMBDAH2
87  self%lcnst%TurbParams(6) = 3000.0_kind_real !ZKMENV
88  self%lcnst%TurbParams(7) = 3000.0_kind_real !ZKHENV
89  self%lcnst%TurbParams(8) = 0.1_kind_real !MINTHICK
90  self%lcnst%TurbParams(9) = 0.0030_kind_real !MINSHEAR
91  self%lcnst%TurbParams(10) = 2.5101471e-8_kind_real !C_B
92  self%lcnst%TurbParams(11) = 1500.0_kind_real !LAMBDA_B
93  self%lcnst%TurbParams(12) = 500.0_kind_real !AKHMMAX
94  self%lcnst%TurbParams(13) = 1.0_kind_real !PRANDTLSFC
95  self%lcnst%TurbParams(14) = 0.75_kind_real !PRANDTLRAD
96  self%lcnst%TurbParams(15) = 0.50_kind_real !BETA_RAD
97  self%lcnst%TurbParams(16) = 0.25_kind_real !BETA_SURF
98  self%lcnst%TurbParams(17) = 0.85_kind_real !KHRADFAC
99  self%lcnst%TurbParams(18) = 0.45_kind_real !KHSFCFAC
100  self%lcnst%TurbParams(19) = 20.0_kind_real !TPFAC_SURF
101  self%lcnst%TurbParams(20) = 1.5e-3_kind_real !ENTRATE_SURF
102  self%lcnst%TurbParams(21) = 0.5_kind_real !PCEFF_SURF
103  self%lcnst%TurbParams(22) = -999.0_kind_real !LOUIS_MEMORY
104  self%lcnst%TurbParamsI(1) = count(pref < 50000.0_kind_real) !KPBLMIN
105  self%lcnst%TurbParamsI(2) = 1 !LOCK_ON
106  self%lcnst%TurbParamsI(3) = 1 !PBLHT_OPTION
107  self%lcnst%TurbParamsI(4) = 0 !RADLW_DEP
108 
109  deallocate(pref)
110 
111 endsubroutine create
112 
113 ! ------------------------------------------------------------------------------
114 
115 subroutine init_nl(self,pert,traj)
117  implicit none
118 
119  class(fv3jedi_lm_turbulence_type), intent(inout) :: self
120  type(fv3jedi_lm_pert), intent(inout) :: pert
121  type(fv3jedi_lm_traj), intent(in) :: traj
122 
123 endsubroutine init_nl
124 
125 ! ------------------------------------------------------------------------------
126 
127 subroutine init_tl(self,pert,traj)
129  implicit none
130 
131  class(fv3jedi_lm_turbulence_type), intent(inout) :: self
132  type(fv3jedi_lm_pert), intent(inout) :: pert
133  type(fv3jedi_lm_traj), intent(in) :: traj
134 
135 endsubroutine init_tl
136 
137 ! ------------------------------------------------------------------------------
138 
139 subroutine init_ad(self,pert,traj)
141  implicit none
142 
143  class(fv3jedi_lm_turbulence_type), intent(inout) :: self
144  type(fv3jedi_lm_pert), intent(inout) :: pert
145  type(fv3jedi_lm_traj), intent(in) :: traj
146 
147 endsubroutine init_ad
148 
149 ! ------------------------------------------------------------------------------
150 
151 subroutine step_nl(self,conf,traj)
153  implicit none
154 
155  class(fv3jedi_lm_turbulence_type), target, intent(inout) :: self
156  type(fv3jedi_lm_traj), target, intent(inout) :: traj
157  type(fv3jedi_lm_conf), intent(in) :: conf
158 
159  type(local_traj_turbulence), pointer :: ltraj
160  real(kind_real), pointer, dimension(:,:,:) :: p_u
161  real(kind_real), pointer, dimension(:,:,:) :: p_v
162  real(kind_real), pointer, dimension(:,:,:) :: p_t
163  real(kind_real), pointer, dimension(:,:,:) :: p_delp
164  real(kind_real), pointer, dimension(:,:,:) :: p_qv
165  real(kind_real), pointer, dimension(:,:,:) :: p_qi
166  real(kind_real), pointer, dimension(:,:,:) :: p_ql
167  real(kind_real), pointer, dimension(:,:,:) :: p_o3
168 
169  !Pointers with ind starting at 1
170  p_u(1:conf%im,1:conf%jm,1:conf%lm) => traj%u
171  p_v(1:conf%im,1:conf%jm,1:conf%lm) => traj%v
172  p_t(1:conf%im,1:conf%jm,1:conf%lm) => traj%t
173  p_delp(1:conf%im,1:conf%jm,1:conf%lm) => traj%delp
174  p_qv(1:conf%im,1:conf%jm,1:conf%lm) => traj%qv
175  p_qi(1:conf%im,1:conf%jm,1:conf%lm) => traj%qi
176  p_ql(1:conf%im,1:conf%jm,1:conf%lm) => traj%ql
177  p_o3(1:conf%im,1:conf%jm,1:conf%lm) => traj%o3
178 
179  !Convenience pointers
180  if (conf%saveltraj) then
181  ltraj => self%ltraj(conf%n)
182  else
183  ltraj => self%ltraj(1)
184  endif
185 
186  !Set up the local trajectory
187  if (.not. ltraj%set) call set_ltraj(conf,self%lcnst,traj,ltraj)
188 
189  !t2pt
190  p_t = p00**kappa * p_t / ltraj%pk
191 
192  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akv,ltraj%bkv,ltraj%ckv,p_u ,1,1)
193  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akv,ltraj%bkv,ltraj%ckv,p_v ,1,1)
194  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%aks,ltraj%bks,ltraj%cks,p_t ,1,1)
195  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_qv,1,1)
196  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_qi,1,0)
197  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_ql,1,0)
198  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_o3,1,0)
199 
200  !pt2t
201  p_t = ltraj%pk * p_t / p00**kappa
202 
203  !Nullify
204  nullify(p_u)
205  nullify(p_v)
206  nullify(p_t)
207  nullify(p_delp)
208  nullify(p_qv)
209  nullify(p_qi)
210  nullify(p_ql)
211  nullify(p_o3)
212  nullify(ltraj)
213 
214 endsubroutine step_nl
215 
216 ! ------------------------------------------------------------------------------
217 
218 subroutine step_tl(self,conf,traj,pert)
220  implicit none
221 
222  class(fv3jedi_lm_turbulence_type), target, intent(inout) :: self
223  type(fv3jedi_lm_conf), intent(in) :: conf
224  type(fv3jedi_lm_traj), intent(in) :: traj
225  type(fv3jedi_lm_pert), target, intent(inout) :: pert
226 
227  type(local_traj_turbulence), pointer :: ltraj
228  real(kind_real), pointer, dimension(:,:,:) :: p_u
229  real(kind_real), pointer, dimension(:,:,:) :: p_v
230  real(kind_real), pointer, dimension(:,:,:) :: p_t
231  real(kind_real), pointer, dimension(:,:,:) :: p_delp
232  real(kind_real), pointer, dimension(:,:,:) :: p_qv
233  real(kind_real), pointer, dimension(:,:,:) :: p_qi
234  real(kind_real), pointer, dimension(:,:,:) :: p_ql
235  real(kind_real), pointer, dimension(:,:,:) :: p_o3
236 
237  !Pointers with ind starting at 1
238  p_u(1:conf%im,1:conf%jm,1:conf%lm) => pert%u
239  p_v(1:conf%im,1:conf%jm,1:conf%lm) => pert%v
240  p_t(1:conf%im,1:conf%jm,1:conf%lm) => pert%t
241  p_delp(1:conf%im,1:conf%jm,1:conf%lm) => pert%delp
242  p_qv(1:conf%im,1:conf%jm,1:conf%lm) => pert%qv
243  p_qi(1:conf%im,1:conf%jm,1:conf%lm) => pert%qi
244  p_ql(1:conf%im,1:conf%jm,1:conf%lm) => pert%ql
245  p_o3(1:conf%im,1:conf%jm,1:conf%lm) => pert%o3
246 
247  !Convenience pointers
248  if (conf%saveltraj) then
249  ltraj => self%ltraj(conf%n)
250  else
251  ltraj => self%ltraj(1)
252  endif
253 
254  !Set up the local trajectory
255  if (.not. ltraj%set) call set_ltraj(conf,self%lcnst,traj,ltraj)
256 
257  !t2pt
258  p_t = p00**kappa * p_t / ltraj%pk
259 
260  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akv,ltraj%bkv,ltraj%ckv,p_u ,1,1)
261  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akv,ltraj%bkv,ltraj%ckv,p_v ,1,1)
262  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%aks,ltraj%bks,ltraj%cks,p_t ,1,1)
263  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_qv,1,1)
264  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_qi,1,0)
265  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_ql,1,0)
266  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_o3,1,0)
267 
268  !pt2t
269  p_t = ltraj%pk * p_t / p00**kappa
270 
271  !Nullify
272  nullify(p_u)
273  nullify(p_v)
274  nullify(p_t)
275  nullify(p_delp)
276  nullify(p_qv)
277  nullify(p_qi)
278  nullify(p_ql)
279  nullify(p_o3)
280  nullify(ltraj)
281 
282 endsubroutine step_tl
283 
284 ! ------------------------------------------------------------------------------
285 
286 subroutine step_ad(self,conf,traj,pert)
288  implicit none
289 
290  class(fv3jedi_lm_turbulence_type), target, intent(inout) :: self
291  type(fv3jedi_lm_conf), intent(in) :: conf
292  type(fv3jedi_lm_traj), intent(in) :: traj
293  type(fv3jedi_lm_pert), target, intent(inout) :: pert
294 
295  type(local_traj_turbulence), pointer :: ltraj
296  real(kind_real), pointer, dimension(:,:,:) :: p_u
297  real(kind_real), pointer, dimension(:,:,:) :: p_v
298  real(kind_real), pointer, dimension(:,:,:) :: p_t
299  real(kind_real), pointer, dimension(:,:,:) :: p_delp
300  real(kind_real), pointer, dimension(:,:,:) :: p_qv
301  real(kind_real), pointer, dimension(:,:,:) :: p_qi
302  real(kind_real), pointer, dimension(:,:,:) :: p_ql
303  real(kind_real), pointer, dimension(:,:,:) :: p_o3
304 
305  !Pointers with ind starting at 1
306  p_u(1:conf%im,1:conf%jm,1:conf%lm) => pert%u
307  p_v(1:conf%im,1:conf%jm,1:conf%lm) => pert%v
308  p_t(1:conf%im,1:conf%jm,1:conf%lm) => pert%t
309  p_delp(1:conf%im,1:conf%jm,1:conf%lm) => pert%delp
310  p_qv(1:conf%im,1:conf%jm,1:conf%lm) => pert%qv
311  p_qi(1:conf%im,1:conf%jm,1:conf%lm) => pert%qi
312  p_ql(1:conf%im,1:conf%jm,1:conf%lm) => pert%ql
313  p_o3(1:conf%im,1:conf%jm,1:conf%lm) => pert%o3
314 
315  !Convenience pointers
316  if (conf%saveltraj) then
317  ltraj => self%ltraj(conf%n)
318  else
319  ltraj => self%ltraj(1)
320  endif
321 
322  !Set up the local trajectory
323  if (.not. ltraj%set) call set_ltraj(conf,self%lcnst,traj,ltraj)
324 
325  !t2pt adjoint
326  p_t = ltraj%pk * p_t / p00**kappa
327 
328  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akv,ltraj%bkv,ltraj%ckv,p_u ,2,1)
329  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akv,ltraj%bkv,ltraj%ckv,p_v ,2,1)
330  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%aks,ltraj%bks,ltraj%cks,p_t ,2,1)
331  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_qv,2,1)
332  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_qi,2,0)
333  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_ql,2,0)
334  call vtrisolvepert(conf%im,conf%jm,conf%lm,ltraj%akq,ltraj%bkq,ltraj%ckq,p_o3,2,0)
335 
336  !pt2t adjoint
337  p_t = p00**kappa * p_t / ltraj%pk
338 
339  !Nullify
340  nullify(p_u)
341  nullify(p_v)
342  nullify(p_t)
343  nullify(p_delp)
344  nullify(p_qv)
345  nullify(p_qi)
346  nullify(p_ql)
347  nullify(p_o3)
348  nullify(ltraj)
349 
350 endsubroutine step_ad
351 
352 ! ------------------------------------------------------------------------------
353 
354 subroutine delete(self,conf)
356  implicit none
357  class(fv3jedi_lm_turbulence_type), intent(inout) :: self
358  type(fv3jedi_lm_conf), intent(in) :: conf
359 
360  integer :: n
361 
362  if (conf%saveltraj) then
363  do n = 1,conf%nt
364  call deallocate_ltraj(self%ltraj(n))
365  enddo
366  else
367  call deallocate_ltraj(self%ltraj(1))
368  endif
369 
370  deallocate(self%ltraj)
371 
372 endsubroutine delete
373 
374 ! ------------------------------------------------------------------------------
375 
376 subroutine set_ltraj(conf,lcnst,traj,ltraj)
378  type(fv3jedi_lm_conf), intent(in) :: conf
379  type(local_cnst_turbulence), intent(in) :: lcnst
380  type(fv3jedi_lm_traj), target, intent(in) :: traj
381  type(local_traj_turbulence), intent(inout) :: ltraj
382 
383  integer :: i,j,l,im,jm,lm,isc,iec,jsc,jec
384 
385  real(kind_real), allocatable, dimension(:,:,:) :: PTT1
386  real(kind_real), allocatable, dimension(:,:) :: ZPBL1, CT1
387  real(kind_real), allocatable, dimension(:,:,:) :: EKV, FKV
388  real(kind_real), allocatable, dimension(:,:,:) :: pet, pmt
389  real(kind_real), allocatable, dimension(:,:,:) :: QIT1, QLT1
390  real(kind_real), allocatable, dimension(:,:,:) :: fQi
391 
392  real(kind_real), pointer, dimension(:,:,:) :: p_u
393  real(kind_real), pointer, dimension(:,:,:) :: p_v
394  real(kind_real), pointer, dimension(:,:,:) :: p_t
395  real(kind_real), pointer, dimension(:,:,:) :: p_delp
396  real(kind_real), pointer, dimension(:,:,:) :: p_qv
397  real(kind_real), pointer, dimension(:,:) :: p_frland
398  real(kind_real), pointer, dimension(:,:) :: p_frocean
399  real(kind_real), pointer, dimension(:,:) :: p_varflt
400  real(kind_real), pointer, dimension(:,:) :: p_cm
401  real(kind_real), pointer, dimension(:,:) :: p_cq
402  real(kind_real), pointer, dimension(:,:) :: p_ustar
403  real(kind_real), pointer, dimension(:,:) :: p_bstar
404 
405  im = conf%im
406  jm = conf%jm
407  lm = conf%lm
408 
409  isc = conf%isc
410  iec = conf%iec
411  jsc = conf%jsc
412  jec = conf%jec
413 
414  allocate(ptt1(1:im,1:jm,1:lm))
415  allocate(ekv(1:im,1:jm,1:lm))
416  allocate(fkv(1:im,1:jm,1:lm))
417  allocate(pet(1:im,1:jm,0:lm))
418  allocate(pmt(1:im,1:jm,1:lm))
419  allocate(qit1(1:im,1:jm,1:lm))
420  allocate(qlt1(1:im,1:jm,1:lm))
421  allocate(fqi(1:im,1:jm,1:lm))
422  allocate(zpbl1(1:im,1:jm))
423  allocate(ct1(1:im,1:jm))
424 
425  p_u(1:im,1:jm,1:lm) => traj%u
426  p_v(1:im,1:jm,1:lm) => traj%v
427  p_t(1:im,1:jm,1:lm) => traj%t
428  p_delp(1:im,1:jm,1:lm) => traj%delp
429  p_qv(1:im,1:jm,1:lm) => traj%qv
430  p_frland(1:im,1:jm) => traj%FRLAND
431  p_frocean(1:im,1:jm) => traj%FROCEAN
432  p_varflt(1:im,1:jm) => traj%VARFLT
433  p_cm(1:im,1:jm) => traj%CM
434  p_cq(1:im,1:jm) => traj%CQ
435  p_ustar(1:im,1:jm) => traj%USTAR
436  p_bstar(1:im,1:jm) => traj%BSTAR
437 
438  !Compute pressures from delp
439  call compute_pressures(im,jm,lm,conf%ptop,p_delp,pet,pmt,ltraj%pk)
440 
441  !Use local copied to avoid overwrite
442  zpbl1(1:im,1:jm) = traj%ZPBL(isc:iec,jsc:jec)
443  ct1(1:im,1:jm) = traj%CT (isc:iec,jsc:jec)
444 
445  ptt1 = p00**kappa * p_t / ltraj%pk
446 
447  !Calculate total cloud ice and liquid trajectory
448  if (conf%do_phy_mst == 0) then
449 
450  qit1(1:im,1:jm,:) = traj%QI(isc:iec,jsc:jec,:)
451  qlt1(1:im,1:jm,:) = traj%QL(isc:iec,jsc:jec,:)
452 
453  else
454 
455  do l = 1,lm
456  do j = 1,jm
457  do i = 1,im
458  call icefraction( p_t(i,j,l), fqi(i,j,l) )
459  enddo
460  enddo
461  enddo
462 
463  qit1(1:im,1:jm,:) = (traj%QLS(isc:iec,jsc:jec,:) + traj%QCN(isc:iec,jsc:jec,:)) * fqi(1:im,1:jm,:)
464  qlt1(1:im,1:jm,:) = (traj%QLS(isc:iec,jsc:jec,:) + traj%QCN(isc:iec,jsc:jec,:)) * (1-fqi(1:im,1:jm,:))
465 
466  endif
467 
468  !Initialize tri-diagonal arrays
469  ltraj%AKV = 0.0
470  ltraj%BKV = 0.0
471  ltraj%CKV = 0.0
472  ltraj%AKS = 0.0
473  ltraj%BKS = 0.0
474  ltraj%CKS = 0.0
475  ltraj%AKQ = 0.0
476  ltraj%BKQ = 0.0
477  ltraj%CKQ = 0.0
478 
479  !Call boundary layer routine. These routines will return the lower (AK*), main (BK*),
480  !and upper (CK*) diagonals.
481 
482  call bl_driver( im , &
483  jm , &
484  lm , &
485  conf%DT , &
486  p_u , &
487  p_v , &
488  ptt1 , &
489  p_qv , &
490  pet , &
491  qit1 , &
492  qlt1 , &
493  p_frland , &
494  p_frocean , &
495  p_varflt , &
496  zpbl1 , &
497  p_cm , &
498  ct1 , &
499  p_cq , &
500  lcnst%TURBPARAMS , &
501  lcnst%TURBPARAMSI , &
502  p_ustar , &
503  p_bstar , &
504  ltraj%AKS, ltraj%BKS, ltraj%CKS, &
505  ltraj%AKQ, ltraj%BKQ, ltraj%CKQ, &
506  ltraj%AKV, ltraj%BKV, ltraj%CKV, &
507  ekv, fkv )
508 
509  !Solver part 1, perform LU decomposition
510  call vtrilupert(im,jm,lm,ltraj%AKV,ltraj%BKV,ltraj%CKV)
511  call vtrilupert(im,jm,lm,ltraj%AKS,ltraj%BKS,ltraj%CKS)
512  call vtrilupert(im,jm,lm,ltraj%AKQ,ltraj%BKQ,ltraj%CKQ)
513 
514  deallocate(ptt1 )
515  deallocate(zpbl1)
516  deallocate(ct1 )
517  deallocate(ekv )
518  deallocate(fkv )
519  deallocate(pet )
520  deallocate(pmt )
521  deallocate(qit1 )
522  deallocate(qlt1 )
523  deallocate(fqi )
524 
525  nullify(p_u)
526  nullify(p_v)
527  nullify(p_t)
528  nullify(p_delp)
529  nullify(p_qv)
530  nullify(p_frland)
531  nullify(p_frocean)
532  nullify(p_varflt)
533  nullify(p_cm)
534  nullify(p_cq)
535  nullify(p_ustar)
536  nullify(p_bstar)
537 
538  if (conf%saveltraj) ltraj%set = .true.
539 
540 endsubroutine set_ltraj
541 
542 ! ------------------------------------------------------------------------------
543 
544 subroutine allocate_ltraj(im,jm,lm,ltraj)
546  integer, intent(in) :: im,jm,lm
547  type(local_traj_turbulence), intent(inout) :: ltraj
548 
549  allocate(ltraj%pk(im,jm,lm))
550  allocate(ltraj%akv(im,jm,lm))
551  allocate(ltraj%bkv(im,jm,lm))
552  allocate(ltraj%ckv(im,jm,lm))
553  allocate(ltraj%aks(im,jm,lm))
554  allocate(ltraj%bks(im,jm,lm))
555  allocate(ltraj%cks(im,jm,lm))
556  allocate(ltraj%akq(im,jm,lm))
557  allocate(ltraj%bkq(im,jm,lm))
558  allocate(ltraj%ckq(im,jm,lm))
559 
560 endsubroutine allocate_ltraj
561 
562 ! ------------------------------------------------------------------------------
563 
564 subroutine deallocate_ltraj(ltraj)
566  type(local_traj_turbulence), intent(inout) :: ltraj
567 
568  deallocate(ltraj%pk )
569  deallocate(ltraj%akv)
570  deallocate(ltraj%bkv)
571  deallocate(ltraj%ckv)
572  deallocate(ltraj%aks)
573  deallocate(ltraj%bks)
574  deallocate(ltraj%cks)
575  deallocate(ltraj%akq)
576  deallocate(ltraj%bkq)
577  deallocate(ltraj%ckq)
578 
579 endsubroutine deallocate_ltraj
580 
581 ! ------------------------------------------------------------------------------
582 
583 subroutine vtrilupert(im,jm,lm,a,b,c)
585  !perform lu decomposition of tridiagonal system
586 
587  implicit none
588 
589  integer, intent(in ) :: im, jm, lm
590  real(kind_real), dimension(im,jm,lm), intent(in ) :: c
591  real(kind_real), dimension(im,jm,lm), intent(inout) :: a, b
592 
593  integer :: l
594 
595  b(:,:,1) = 1. / b(:,:,1)
596  do l = 2,lm
597  a(:,:,l) = a(:,:,l) * b(:,:,l-1)
598  b(:,:,l) = 1. / ( b(:,:,l) - c(:,:,l-1) * a(:,:,l) )
599  end do
600 
601 end subroutine vtrilupert
602 
603 ! ------------------------------------------------------------------------------
604 
605 subroutine vtrisolvepert(im,jm,lm,a,b,c,y,phase,ygswitch)
607  !solve lu decomposed tridiagonal system
608 
609  implicit none
610 
611  !arguments
612  integer, intent(in ) :: im, jm, lm, phase, ygswitch
613  real(kind_real), dimension(im,jm,lm), intent(in ) :: a, b, c
614  real(kind_real), dimension(im,jm,lm), intent(inout) :: y
615 
616  !locals
617  integer :: l
618 
619  if (phase == 1) then !TLM
620 
621  !solve (lu)ynew = yold
622 
623  !sweep down, modifying rhs with multiplier a
624  do l=2,lm
625  y(:,:,l) = y(:,:,l) - a(:,:,l)*y(:,:,l-1)
626  enddo
627 
628  !sweep up, solving for updated value. note b has the inverse of the main diagonal
629  if (ygswitch == 1) then !winds, temperature and q
630  y(:,:,lm) = y(:,:,lm)*b(:,:,lm) !ygprime = 0
631  else !tracers
632  y(:,:,lm) = y(:,:,lm)*b(:,:,lm-1)/(b(:,:,lm-1) - a(:,:,lm)*(1.0+c(:,:,lm-1)*b(:,:,lm-1) ))
633  endif
634 
635  do l=lm-1,1,-1
636  y(:,:,l) = b(:,:,l)*(y(:,:,l)-c(:,:,l)*y(:,:,l+1) )
637  enddo
638 
639  elseif (phase == 2) then !ADM
640 
641  if (ygswitch == 1) then !can solve (lu)'ynew = u'l'ynew = yold
642 
643  !sweep down but with u' instead of l
644  y(:,:,1) = y(:,:,1)*b(:,:,1)
645  do l=2,lm
646  y(:,:,l) = b(:,:,l) * (y(:,:,l) - c(:,:,l-1)*y(:,:,l-1))
647  enddo
648 
649  !sweep up but with l' instead of u
650  do l=lm-1,1,-1
651  y(:,:,l) = y(:,:,l) - a(:,:,l+1)*y(:,:,l+1)
652  enddo
653 
654  else !change for surface means transpose doesnt work so use line-by-line adjoint
655 
656  !adjoint of sweep up, solving for updated value. note b has the inverse of the main diagonal
657  do l=1,lm-1,1
658  y(:,:,l+1) = y(:,:,l+1) - c(:,:,l)*b(:,:,l)*y(:,:,l)
659  y(:,:,l) = b(:,:,l)*y(:,:,l)
660  end do
661 
662  !adjoint of surface fix
663  y(:,:,lm) = b(:,:,lm-1)*y(:,:,lm)/(b(:,:,lm-1)-a(:,:,lm)*(c(:,:,lm-1)*b(:,:,lm-1)+1.0))
664 
665  !adjoint of sweep down, modifying rhs with multiplier a
666  do l=lm,2,-1
667  y(:,:,l-1) = y(:,:,l-1) - a(:,:,l)*y(:,:,l)
668  end do
669 
670  endif
671 
672  endif
673 
674 end subroutine vtrisolvepert
675 
676 ! ------------------------------------------------------------------------------
677 
678 end module fv3jedi_lm_turbulence_mod
subroutine allocate_ltraj(im, jm, lm, ltraj)
subroutine vtrilupert(im, jm, lm, a, b, c)
integer, parameter, public set
real(kind=kind_real), parameter, public p00
Compute ice fraction from temperature.
subroutine step_nl(self, conf, traj)
type(cp_iter_controls_type), target, public cp_iter_controls
Definition: conf.py:1
subroutine, public bl_driver(IM, JM, LM, DT, U, V, TH, Q, P, QIT, QLT, FRLAND, FROCEAN, VARFLT, ZPBL, CM, CT, CQ, TURBPARAMS, TURBPARAMSI, USTAR, BSTAR, AKS, BKS, CKS, AKQ, BKQ, CKQ, AKV, BKV, CKV, EKV, FKV)
Definition: bldriver.F90:28
subroutine, public cp_mod_ini(cp_mod_index)
subroutine, public cp_mod_mid
subroutine step_tl(self, conf, traj, pert)
subroutine init_tl(self, pert, traj)
subroutine, public cp_mod_end
subroutine set_ltraj(conf, lcnst, traj, ltraj)
subroutine, public delete(self)
subroutine init_ad(self, pert, traj)
type(cp_iter_type), dimension(:), allocatable, target, public cp_iter
subroutine, public finalize_cp_iter
subroutine step_ad(self, conf, traj, pert)
subroutine init_nl(self, pert, traj)
subroutine, public initialize_cp_iter
subroutine, public create(self, geom, vars)
real(kind=kind_real), parameter, public kappa
Constants for the FV3 model.
subroutine vtrisolvepert(im, jm, lm, a, b, c, y, phase, ygswitch)