Line data Source code
1 : module geopotential_temp
2 :
3 : !---------------------------------------------------------------------------
4 : ! Compute geopotential from temperature
5 : !
6 : ! The hydrostatic matrix elements must be consistent with the dynamics
7 : ! algorithm.
8 : ! The diagonal element is the itegration weight from interface kint minus
9 : ! lyr_step to midpoint kint.
10 : ! The offdiagonal element is the weight between interfaces.
11 : !
12 : ! Author: B.Boville, Feb 2001 from earlier code by Boville and S.J. Lin
13 : !---------------------------------------------------------------------------
14 :
15 : use ccpp_kinds, only: kind_phys
16 : use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
17 :
18 : implicit none
19 : private
20 : save
21 :
22 : public geopotential_temp_run
23 :
24 : CONTAINS
25 : !===========================================================================
26 :
27 : !> \section arg_table_geopotential_temp_run Argument Table
28 : !! \htmlinclude geopotential_temp_run.html
29 23907312 : subroutine geopotential_temp_run(pver, lagrang, layer_surf, layer_toa, &
30 23907312 : interface_surf, interface_toa, ncnst, piln, pint, pmid, pdel, rpdel, &
31 47814624 : temp, qv, carr, cprops, rair, gravit, zvir, zi, zm, ncol, &
32 0 : errflg, errmsg)
33 :
34 : !-----------------------------------------------------------------------
35 : !
36 : ! Purpose:
37 : ! Compute the geopotential height (above the surface) at the midpoints and
38 : ! interfaces using the input temperatures and pressures.
39 : !
40 : !-----------------------------------------------------------------------
41 :
42 : !------------------------------Arguments--------------------------------
43 : !
44 : ! Input arguments
45 : !
46 : integer, intent(in) :: pver
47 : logical, intent(in) :: lagrang ! lagrangian vertical coordinate?
48 : integer, intent(in) :: layer_surf
49 : integer, intent(in) :: layer_toa
50 : integer, intent(in) :: interface_surf
51 : integer, intent(in) :: interface_toa
52 : integer, intent(in) :: ncnst ! Number of constituents
53 : integer, intent(in) :: ncol ! Number of columns
54 :
55 : ! piln: Log interface pressures
56 : real(kind_phys), intent(in) :: piln(:,:) ! (ncol,pverp)
57 : ! pint: Interface pressures
58 : real(kind_phys), intent(in) :: pint(:,:) ! (ncol,pverp)
59 : ! pmid: Midpoint pressures
60 : real(kind_phys), intent(in) :: pmid(:,:) ! (ncol,pver)
61 : ! pdel: layer thickness
62 : real(kind_phys), intent(in) :: pdel(:,:) ! (ncol,pver)
63 : ! rpdel: inverse of layer thickness
64 : real(kind_phys), intent(in) :: rpdel(:,:) ! (ncol,pver)
65 : ! temp: temperature
66 : real(kind_phys), intent(in) :: temp(:,:) ! (ncol,pver)
67 : ! qv: specific humidity
68 : real(kind_phys), intent(in) :: qv(:,:) ! (ncol,pver)
69 : ! carr: ccpp constituents
70 : real(kind_phys), intent(in) :: carr(:,:,:) ! (ncol, pver, num_const)
71 : ! cprops: ccpp constituent properties
72 : type(ccpp_constituent_prop_ptr_t), intent(in) :: cprops(:) ! (num_const)
73 : ! rair: Gas constant for dry air
74 : real(kind_phys), intent(in) :: rair(:,:) ! (ncol,pver)
75 : ! gravit: Acceleration of gravity
76 : real(kind_phys), intent(in) :: gravit
77 : ! zvir: rh2o/rair - 1
78 : real(kind_phys), intent(in) :: zvir(:,:) ! (ncol,pver)
79 :
80 : ! Output arguments
81 :
82 : ! zi: Height above surface at interfaces
83 : real(kind_phys), intent(out) :: zi(:,:) ! (ncol,pverp)
84 : ! zm: Geopotential height at mid level
85 : real(kind_phys), intent(out) :: zm(:,:) ! (ncol,pver)
86 : integer, intent(out) :: errflg
87 : character(len=512), intent(out) :: errmsg
88 : !
89 : !---------------------------Local variables-----------------------------
90 : !
91 : logical :: active_flag ! Flag for whether consitutent is thermodynamically active
92 : logical :: water_flag ! Flag for whether constituent is a water species
93 : integer :: icol ! Horizontal index
94 : integer :: klyr ! Vertical layer index
95 : integer :: kint ! Vertical interface index
96 : integer :: cidx ! CCPP constituent index
97 : integer :: lyr_step ! increment to move up a level
98 : integer :: int_step ! increment to move up a level
99 47814624 : real(kind_phys) :: hkk(ncol) ! diagonal element of hydrostatic matrix
100 47814624 : real(kind_phys) :: hkl(ncol) ! off-diagonal element
101 47814624 : real(kind_phys) :: rog(ncol,pver) ! Rair / gravit
102 : real(kind_phys) :: tv ! virtual temperature
103 : real(kind_phys) :: tvfac ! Tv / T
104 47814624 : real(kind_phys) :: qfac(ncol,pver) ! factor to convert from wet to dry mixing ratio
105 23907312 : real(kind_phys) :: sum_dry_mixing_ratio(ncol,pver) ! sum of dry water mixing ratios
106 : !
107 : !-----------------------------------------------------------------------
108 : !
109 :
110 23907312 : errmsg = ''
111 23907312 : errflg = 0
112 :
113 10403016624 : rog(:ncol,:) = rair(:ncol,:) / gravit
114 :
115 : ! Determine order of vertical dimension loop
116 23907312 : if (layer_surf > layer_toa) then
117 : lyr_step = -1
118 : else
119 0 : lyr_step = 1
120 : end if
121 23907312 : if (interface_surf > interface_toa) then
122 : int_step = -1
123 : else
124 0 : int_step = 1
125 : end if
126 :
127 : ! The surface height is zero by definition.
128 399196512 : zi(:ncol, interface_surf) = 0.0_kind_phys
129 :
130 : ! For the computation of generalized virtual temperature (equation 16
131 : ! in Lauritzen et al. (2018); https://doi.org/10.1029/2017MS001257)
132 : !
133 : ! Compute factor for converting wet to dry mixing ratio (eq.7)
134 : !
135 10403016624 : qfac = 1._kind_phys
136 23907312 : active_flag = .false.
137 23907312 : water_flag = .false.
138 95629248 : do cidx = 1,ncnst
139 : ! Check if constituent is thermodynamically active:
140 71721936 : call cprops(cidx)%is_thermo_active(active_flag)
141 95629248 : if (active_flag) then
142 : ! Check if consitutent is also a water species:
143 71721936 : call cprops(cidx)%is_water_species(water_flag)
144 71721936 : if (water_flag) then
145 1936492272 : do klyr = layer_surf, layer_toa, lyr_step
146 31209049872 : do icol = 1, ncol
147 : ! Add thermodynamically active water species
148 : ! to mixing ratio factor:
149 58545115200 : qfac(icol, klyr) = qfac(icol, klyr) - &
150 89682443136 : carr(icol, klyr, cidx)
151 : end do
152 : end do
153 : end if
154 : end if
155 : end do
156 10403016624 : qfac = 1._kind_phys/qfac
157 :
158 : ! Compute sum of thermodynamically active water
159 : ! constituent mixing ratios with respect to dry air:
160 10403016624 : sum_dry_mixing_ratio = 1._kind_phys
161 23907312 : active_flag = .false.
162 23907312 : water_flag = .false.
163 95629248 : do cidx = 1,ncnst
164 : !Check if constituent is thermodynamically active:
165 71721936 : call cprops(cidx)%is_thermo_active(active_flag)
166 95629248 : if (active_flag) then
167 : !Check if consitutent is also a water species:
168 71721936 : call cprops(cidx)%is_water_species(water_flag)
169 71721936 : if (water_flag) then
170 1936492272 : do klyr = layer_surf, layer_toa, lyr_step
171 31209049872 : do icol = 1, ncol
172 58545115200 : sum_dry_mixing_ratio(icol, klyr) = &
173 : sum_dry_mixing_ratio(icol, klyr) + &
174 89682443136 : carr(icol, klyr, cidx) * qfac(icol,klyr)
175 : end do
176 : end do
177 : end if
178 : end if
179 : end do
180 10403016624 : sum_dry_mixing_ratio(:,:) = 1._kind_phys/sum_dry_mixing_ratio(:,:)
181 :
182 : ! Compute zi, zm from bottom up.
183 : ! Note, zi(i,k) is the interface above zm(i,k) when ordered top to bot
184 :
185 23907312 : klyr = layer_surf
186 23907312 : do kint = interface_surf + int_step, interface_toa, int_step
187 :
188 : ! First set hydrostatic elements consistent with dynamics
189 621590112 : if (lagrang) then
190 0 : do icol = 1, ncol
191 0 : hkl(icol) = piln(icol, kint-int_step) - piln(icol,kint)
192 : hkk(icol) = 1._kind_phys - &
193 0 : (pint(icol,kint) * hkl(icol) * rpdel(icol,klyr))
194 : end do
195 : else
196 10379109312 : do icol = 1,ncol
197 9757519200 : hkl(icol) = pdel(icol,klyr) / pmid(icol,klyr)
198 10379109312 : hkk(icol) = 0.5_kind_phys * hkl(icol)
199 : end do
200 : end if
201 :
202 : ! Now compute tv, zm, zi
203 :
204 10379109312 : do icol = 1, ncol
205 19515038400 : tvfac = (1._kind_phys + (zvir(icol,klyr) + 1._kind_phys) * &
206 9757519200 : qv(icol,klyr)*qfac(icol,klyr)) * &
207 19515038400 : sum_dry_mixing_ratio(icol,klyr)
208 9757519200 : tv = temp(icol,klyr) * tvfac
209 :
210 19515038400 : zm(icol,klyr) = zi(icol,kint-int_step) + &
211 9757519200 : (rog(icol,klyr) * tv * hkk(icol))
212 9757519200 : zi(icol,kint) = zi(icol,kint-int_step) + &
213 20136628512 : (rog(icol,klyr) * tv * hkl(icol))
214 : end do
215 621590112 : klyr = klyr + lyr_step
216 : end do
217 :
218 23907312 : end subroutine geopotential_temp_run
219 :
220 : end module geopotential_temp
|
