Line data Source code
1 : module dadadj
2 : !=======================================================================
3 : ! GFDL style dry adiabatic adjustment
4 : !
5 : ! Method:
6 : ! if stratification is unstable, adjustment to the dry adiabatic lapse
7 : ! rate is forced subject to the condition that enthalpy is conserved.
8 : !=======================================================================
9 :
10 : use ccpp_kinds, only: kind_phys
11 :
12 : implicit none
13 : private
14 : save
15 :
16 : public :: dadadj_init ! init routine
17 : public :: dadadj_run ! main routine
18 :
19 : integer :: nlvdry ! number of layers from top of model to apply the adjustment
20 : integer :: niter ! number of iterations for convergence
21 :
22 : CONTAINS
23 :
24 : !> \section arg_table_dadadj_init Argument Table
25 : !! \htmlinclude dadadj_init.html
26 1536 : subroutine dadadj_init(dadadj_nlvdry, dadadj_niter, nz, errmsg, errflg)
27 : !------------------------------------------------
28 : ! Input / output parameters
29 : !------------------------------------------------
30 : integer, intent(in) :: dadadj_nlvdry
31 : integer, intent(in) :: dadadj_niter
32 : integer, intent(in) :: nz
33 : character(len=512), intent(out) :: errmsg
34 : integer, intent(out) :: errflg
35 :
36 1536 : errmsg = ''
37 1536 : errflg = 0
38 :
39 1536 : if (dadadj_nlvdry >= nz .or. dadadj_nlvdry < 0) then
40 0 : errflg = 1
41 0 : write(errmsg,*) 'dadadj_init: dadadj_nlvdry=',dadadj_nlvdry,' but must be less than the number of vertical levels ',&
42 0 : '(',nz,'), and must be a positive integer.'
43 : end if
44 :
45 1536 : nlvdry = dadadj_nlvdry
46 1536 : niter = dadadj_niter
47 :
48 1536 : end subroutine dadadj_init
49 :
50 : !> \section arg_table_dadadj_run Argument Table
51 : !! \htmlinclude dadadj_run.html
52 1495368 : subroutine dadadj_run( &
53 2990736 : ncol, nz, dt, pmid, pint, pdel, state_t, state_q, cappa, cpair, s_tend, &
54 4486104 : q_tend, dadpdf, scheme_name, errmsg, errflg)
55 :
56 : !------------------------------------------------
57 : ! Input / output parameters
58 : !------------------------------------------------
59 : integer, intent(in) :: ncol ! number of atmospheric columns
60 : integer, intent(in) :: nz ! number of atmospheric levels
61 : real(kind_phys), intent(in) :: dt ! physics timestep
62 : real(kind_phys), intent(in) :: pmid(:,:) ! pressure at model levels
63 : real(kind_phys), intent(in) :: pint(:,:) ! pressure at model interfaces
64 : real(kind_phys), intent(in) :: pdel(:,:) ! vertical delta-p
65 : real(kind_phys), intent(in) :: cappa(:,:) ! variable Kappa
66 : real(kind_phys), intent(in) :: cpair(:,:) ! heat capacity of air
67 : real(kind_phys), intent(in) :: state_t(:,:) ! temperature (K)
68 : real(kind_phys), intent(in) :: state_q(:,:) ! specific humidity
69 : real(kind_phys), intent(out) :: s_tend(:,:) ! dry air enthalpy tendency
70 : real(kind_phys), intent(out) :: q_tend(:,:) ! specific humidity tendency
71 : real(kind_phys), intent(out) :: dadpdf(:,:) ! PDF of where adjustments happened
72 :
73 : character(len=64), intent(out) :: scheme_name
74 : character(len=512), intent(out) :: errmsg
75 : integer, intent(out) :: errflg
76 :
77 : !------------------------------------------------
78 : ! Local variables
79 : !------------------------------------------------
80 :
81 : integer :: i,k ! longitude, level indices
82 : integer :: jiter ! iteration index
83 1495368 : real(kind_phys), allocatable :: c1dad(:) ! intermediate constant
84 1495368 : real(kind_phys), allocatable :: c2dad(:) ! intermediate constant
85 1495368 : real(kind_phys), allocatable :: c3dad(:) ! intermediate constant
86 1495368 : real(kind_phys), allocatable :: c4dad(:) ! intermediate constant
87 : real(kind_phys) :: gammad ! dry adiabatic lapse rate (deg/Pa)
88 : real(kind_phys) :: zeps ! convergence criterion (deg/Pa)
89 : real(kind_phys) :: rdenom ! reciprocal of denominator of expression
90 : real(kind_phys) :: dtdp ! delta-t/delta-p
91 : real(kind_phys) :: zepsdp ! zeps*delta-p
92 : real(kind_phys) :: zgamma ! intermediate constant
93 : real(kind_phys) :: qave ! mean q between levels
94 : real(kind_phys) :: cappaint ! Kappa at level intefaces
95 2990736 : real(kind_phys) :: t(ncol,nz)
96 2990736 : real(kind_phys) :: q(ncol,nz)
97 :
98 : logical :: ilconv ! .TRUE. ==> convergence was attained
99 1495368 : logical :: dodad(ncol) ! .TRUE. ==> do dry adjustment
100 :
101 : !-----------------------------------------------------------------------
102 :
103 1495368 : zeps = 2.0e-5_kind_phys ! set convergence criteria
104 1495368 : errmsg = ''
105 1495368 : errflg = 0
106 1495368 : scheme_name = 'DADADJ'
107 :
108 4486104 : allocate(c1dad(nlvdry), stat=errflg)
109 1495368 : if (errflg /= 0) then
110 0 : errmsg = trim(scheme_name)//': Allocate of c1dad(nlvdry) failed'
111 0 : return
112 : end if
113 4486104 : allocate(c2dad(nlvdry), stat=errflg)
114 1495368 : if (errflg /= 0) then
115 0 : errmsg = trim(scheme_name)//': Allocate of c2dad(nlvdry) failed'
116 0 : return
117 : end if
118 4486104 : allocate(c3dad(nlvdry), stat=errflg)
119 1495368 : if (errflg /= 0) then
120 0 : errmsg = trim(scheme_name)//': Allocate of c3dad(nlvdry) failed'
121 0 : return
122 : end if
123 4486104 : allocate(c4dad(nlvdry), stat=errflg)
124 1495368 : if (errflg /= 0) then
125 0 : errmsg = trim(scheme_name)//': Allocate of c4dad(nlvdry) failed'
126 0 : return
127 : end if
128 :
129 2323627992 : t = state_t
130 2323627992 : q = state_q
131 :
132 : ! Find gridpoints with unstable stratification
133 :
134 24969168 : do i = 1, ncol
135 23473800 : cappaint = 0.5_kind_phys*(cappa(i,2) + cappa(i,1))
136 23473800 : gammad = cappaint*0.5_kind_phys*(t(i,2) + t(i,1))/pint(i,2)
137 23473800 : dtdp = (t(i,2) - t(i,1))/(pmid(i,2) - pmid(i,1))
138 24969168 : dodad(i) = (dtdp + zeps) > gammad
139 : end do
140 :
141 2323627992 : dadpdf(:ncol,:) = 0._kind_phys
142 4486104 : do k= 2, nlvdry
143 51433704 : do i = 1, ncol
144 46947600 : cappaint = 0.5_kind_phys*(cappa(i,k+1) + cappa(i,k))
145 46947600 : gammad = cappaint*0.5_kind_phys*(t(i,k+1) + t(i,k))/pint(i,k+1)
146 46947600 : dtdp = (t(i,k+1) - t(i,k))/(pmid(i,k+1) - pmid(i,k))
147 46947600 : dodad(i) = dodad(i) .or. (dtdp + zeps) > gammad
148 49938336 : if ((dtdp + zeps) > gammad) then
149 19402 : dadpdf(i,k) = 1._kind_phys
150 : end if
151 : end do
152 : end do
153 :
154 : ! Make a dry adiabatic adjustment
155 : ! Note: nlvdry ****MUST**** be < pver
156 :
157 24969168 : COL: do i = 1, ncol
158 :
159 24969168 : if (dodad(i)) then
160 :
161 89577 : zeps = 2.0e-5_kind_phys
162 :
163 358308 : do k = 1, nlvdry
164 268731 : cappaint = 0.5_kind_phys*(cappa(i,k+1) + cappa(i,k))
165 268731 : c1dad(k) = cappaint*0.5_kind_phys*(pmid(i,k+1)-pmid(i,k))/pint(i,k+1)
166 268731 : c2dad(k) = (1._kind_phys - c1dad(k))/(1._kind_phys + c1dad(k))
167 268731 : rdenom = 1._kind_phys/(pdel(i,k)*c2dad(k) + pdel(i,k+1))
168 268731 : c3dad(k) = rdenom*pdel(i,k)
169 358308 : c4dad(k) = rdenom*pdel(i,k+1)
170 : end do
171 :
172 : ilconv = .false.
173 :
174 89577 : DBLZEP: do while (.not. ilconv)
175 :
176 182793 : do jiter = 1, niter
177 : ilconv = .true.
178 :
179 731172 : do k = 1, nlvdry
180 548379 : zepsdp = zeps*(pmid(i,k+1) - pmid(i,k))
181 548379 : zgamma = c1dad(k)*(t(i,k) + t(i,k+1))
182 :
183 731172 : if ((t(i,k+1)-t(i,k)) >= (zgamma+zepsdp)) then
184 96672 : ilconv = .false.
185 96672 : t(i,k+1) = t(i,k)*c3dad(k) + t(i,k+1)*c4dad(k)
186 96672 : t(i,k) = c2dad(k)*t(i,k+1)
187 96672 : qave = (pdel(i,k+1)*q(i,k+1) + pdel(i,k)*q(i,k))/(pdel(i,k+1)+ pdel(i,k))
188 96672 : q(i,k+1) = qave
189 96672 : q(i,k) = qave
190 : end if
191 :
192 : end do
193 :
194 182793 : if (ilconv) cycle COL ! convergence => next longitude
195 :
196 : end do
197 :
198 0 : zeps = zeps + zeps
199 0 : if (zeps > 1.e-4_kind_phys) then
200 0 : errflg = i
201 0 : write(errmsg,*) trim(scheme_name)//': Convergence failure at column ',i,' zeps > 1.e-4 '// &
202 0 : '(errflg set to failing column index)'
203 0 : return ! error return
204 : end if
205 : end do DBLZEP
206 :
207 : end if
208 :
209 : end do COL
210 :
211 2323627992 : s_tend = (t - state_t)/dt*cpair
212 2323627992 : q_tend = (q - state_q)/dt
213 :
214 1495368 : deallocate(c1dad, c2dad, c3dad, c4dad)
215 :
216 1495368 : end subroutine dadadj_run
217 :
218 : end module dadadj
|
