Line data Source code
1 : module dadadj
2 : !-----------------------------------------------------------------------
3 : !
4 : ! Purpose:
5 : ! GFDL style dry adiabatic adjustment
6 : !
7 : ! Method:
8 : ! if stratification is unstable, adjustment to the dry adiabatic lapse
9 : ! rate is forced subject to the condition that enthalpy is conserved.
10 : !
11 : ! Author: J.Hack
12 : !
13 : !-----------------------------------------------------------------------
14 :
15 : use shr_kind_mod, only: r8 => shr_kind_r8
16 :
17 : implicit none
18 : private
19 : save
20 :
21 : public :: &
22 : dadadj_initial, &
23 : dadadj_calc
24 :
25 : integer :: nlvdry ! number of layers from top of model to apply the adjustment
26 : integer :: niter ! number of iterations for convergence
27 :
28 : !===============================================================================
29 : contains
30 : !===============================================================================
31 :
32 1536 : subroutine dadadj_initial(nlvdry_in, niter_in)
33 :
34 : integer, intent(in) :: nlvdry_in
35 : integer, intent(in) :: niter_in
36 :
37 1536 : nlvdry = nlvdry_in
38 1536 : niter = niter_in
39 :
40 1536 : end subroutine dadadj_initial
41 :
42 : !===============================================================================
43 :
44 65016 : subroutine dadadj_calc( &
45 65016 : ncol, pmid, pint, pdel, cappav, t, &
46 65016 : q, dadpdf, icol_err)
47 :
48 : ! Arguments
49 :
50 : integer, intent(in) :: ncol ! number of atmospheric columns
51 :
52 : real(r8), intent(in) :: pmid(:,:) ! pressure at model levels
53 : real(r8), intent(in) :: pint(:,:) ! pressure at model interfaces
54 : real(r8), intent(in) :: pdel(:,:) ! vertical delta-p
55 : real(r8), intent(in) :: cappav(:,:) ! variable Kappa
56 :
57 : real(r8), intent(inout) :: t(:,:) ! temperature (K)
58 : real(r8), intent(inout) :: q(:,:) ! specific humidity
59 :
60 : real(r8), intent(out) :: dadpdf(:,:) ! PDF of where adjustments happened
61 :
62 : integer, intent(out) :: icol_err ! index of column in which error occurred
63 :
64 : !---------------------------Local workspace-----------------------------
65 :
66 : integer :: i,k ! longitude, level indices
67 : integer :: jiter ! iteration index
68 :
69 65016 : real(r8), allocatable :: c1dad(:) ! intermediate constant
70 65016 : real(r8), allocatable :: c2dad(:) ! intermediate constant
71 65016 : real(r8), allocatable :: c3dad(:) ! intermediate constant
72 65016 : real(r8), allocatable :: c4dad(:) ! intermediate constant
73 : real(r8) :: gammad ! dry adiabatic lapse rate (deg/Pa)
74 : real(r8) :: zeps ! convergence criterion (deg/Pa)
75 : real(r8) :: rdenom ! reciprocal of denominator of expression
76 : real(r8) :: dtdp ! delta-t/delta-p
77 : real(r8) :: zepsdp ! zeps*delta-p
78 : real(r8) :: zgamma ! intermediate constant
79 : real(r8) :: qave ! mean q between levels
80 : real(r8) :: cappa ! Kappa at level intefaces
81 :
82 : logical :: ilconv ! .TRUE. ==> convergence was attained
83 130032 : logical :: dodad(ncol) ! .TRUE. ==> do dry adjustment
84 :
85 : !-----------------------------------------------------------------------
86 :
87 65016 : icol_err = 0
88 65016 : zeps = 2.0e-5_r8 ! set convergence criteria
89 :
90 390096 : allocate(c1dad(nlvdry), c2dad(nlvdry), c3dad(nlvdry), c4dad(nlvdry))
91 :
92 : ! Find gridpoints with unstable stratification
93 :
94 1085616 : do i = 1, ncol
95 1020600 : cappa = 0.5_r8*(cappav(i,2) + cappav(i,1))
96 1020600 : gammad = cappa*0.5_r8*(t(i,2) + t(i,1))/pint(i,2)
97 1020600 : dtdp = (t(i,2) - t(i,1))/(pmid(i,2) - pmid(i,1))
98 1085616 : dodad(i) = (dtdp + zeps) .gt. gammad
99 : end do
100 :
101 101027304 : dadpdf(:ncol,:) = 0._r8
102 195048 : do k= 2, nlvdry
103 2236248 : do i = 1, ncol
104 2041200 : cappa = 0.5_r8*(cappav(i,k+1) + cappav(i,k))
105 2041200 : gammad = cappa*0.5_r8*(t(i,k+1) + t(i,k))/pint(i,k+1)
106 2041200 : dtdp = (t(i,k+1) - t(i,k))/(pmid(i,k+1) - pmid(i,k))
107 2041200 : dodad(i) = dodad(i) .or. (dtdp + zeps).gt.gammad
108 2171232 : if ((dtdp + zeps).gt.gammad) then
109 1097 : dadpdf(i,k) = 1._r8
110 : end if
111 : end do
112 : end do
113 :
114 : ! Make a dry adiabatic adjustment
115 : ! Note: nlvdry ****MUST**** be < pver
116 :
117 1085616 : COL: do i = 1, ncol
118 :
119 1085616 : if (dodad(i)) then
120 :
121 9488 : zeps = 2.0e-5_r8
122 :
123 37952 : do k = 1, nlvdry
124 28464 : c1dad(k) = cappa*0.5_r8*(pmid(i,k+1)-pmid(i,k))/pint(i,k+1)
125 28464 : c2dad(k) = (1._r8 - c1dad(k))/(1._r8 + c1dad(k))
126 28464 : rdenom = 1._r8/(pdel(i,k)*c2dad(k) + pdel(i,k+1))
127 28464 : c3dad(k) = rdenom*pdel(i,k)
128 37952 : c4dad(k) = rdenom*pdel(i,k+1)
129 : end do
130 :
131 : 50 continue
132 :
133 19288 : do jiter = 1, niter
134 : ilconv = .true.
135 :
136 77152 : do k = 1, nlvdry
137 57864 : zepsdp = zeps*(pmid(i,k+1) - pmid(i,k))
138 57864 : zgamma = c1dad(k)*(t(i,k) + t(i,k+1))
139 :
140 77152 : if ((t(i,k+1)-t(i,k)) >= (zgamma+zepsdp)) then
141 10111 : ilconv = .false.
142 10111 : t(i,k+1) = t(i,k)*c3dad(k) + t(i,k+1)*c4dad(k)
143 10111 : t(i,k) = c2dad(k)*t(i,k+1)
144 10111 : qave = (pdel(i,k+1)*q(i,k+1) + pdel(i,k)*q(i,k))/(pdel(i,k+1)+ pdel(i,k))
145 10111 : q(i,k+1) = qave
146 10111 : q(i,k) = qave
147 : end if
148 :
149 : end do
150 :
151 19288 : if (ilconv) cycle COL ! convergence => next longitude
152 : end do
153 :
154 : ! Double convergence criterion if no convergence in niter iterations
155 :
156 0 : zeps = zeps + zeps
157 0 : if (zeps > 1.e-4_r8) then
158 0 : icol_err = i
159 0 : return ! error return
160 : else
161 : go to 50
162 : end if
163 :
164 : end if
165 :
166 : end do COL
167 :
168 65016 : deallocate(c1dad, c2dad, c3dad, c4dad)
169 :
170 65016 : end subroutine dadadj_calc
171 :
172 : !===============================================================================
173 :
174 : end module dadadj
|
