CCPP SciDoc v7.0.x  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_h2ophys.F90
1
3
4module module_h2ophys
5 use machine, only : kind_phys
6 implicit none
7
8 public ty_h2ophys
9
10! #########################################################################################
17! #########################################################################################
18 type ty_h2ophys
19 integer :: nlat
20 integer :: nlev
21 integer :: ntime
22 integer :: ncf
23 real(kind_phys), allocatable :: lat(:)
24 real(kind_phys), allocatable :: pres(:)
25 real(kind_phys), allocatable :: ph2o(:)
26 real(kind_phys), allocatable :: time(:)
27 real(kind_phys), allocatable :: data(:,:,:,:)
28 contains
29 procedure, public :: load
30 procedure, public :: setup
31 procedure, public :: update
32 procedure, public :: run
33 end type ty_h2ophys
34
35contains
36 ! #########################################################################################
37 ! Procedure (type-bound) for loading data.
38 ! #########################################################################################
39 function load(this, file, fileID) result (err_message)
40 class(ty_h2ophys), intent(inout) :: this
41 integer, intent(in) :: fileid
42 character(len=*), intent(in) :: file
43 character(len=128) :: err_message
44 integer :: i1, i2, i3, ierr
45 real(kind=4), dimension(:), allocatable :: lat4, pres4, time4, tempin
46
47 ! initialize error message
48 err_message = ""
49
50 ! Get dimensions from data file
51 open(unit=fileid,file=trim(file), form='unformatted', convert='big_endian', iostat=ierr, iomsg=err_message)
52 if (ierr /= 0 ) return
53 read (fileid, iostat=ierr, iomsg=err_message) this%ncf, this%nlat, this%nlev, this%ntime
54 if (ierr /= 0 ) return
55 rewind(fileid)
56
57 allocate (this%lat(this%nlat))
58 allocate (this%pres(this%nlev))
59 allocate (this%ph2o(this%nlev))
60 allocate (this%time(this%ntime+1))
61 allocate (this%data(this%nlat,this%nlev,this%ncf,this%ntime))
62
63 allocate(lat4(this%nlat), pres4(this%nlev), time4(this%ntime+1))
64 read (fileid, iostat=ierr, iomsg=err_message) this%ncf, this%nlat, this%nlev, this%ntime, lat4, pres4, time4
65 if (ierr /= 0 ) return
66
67 ! Store
68 this%pres(:) = pres4(:)
69 this%ph2o(:) = log(100.0*this%pres(:)) ! from mb to ln(Pa)
70 this%lat(:) = lat4(:)
71 this%time(:) = time4(:)
72 deallocate(lat4, pres4, time4)
73
74 allocate(tempin(this%nlat))
75 do i1=1,this%ntime
76 do i2=1,this%ncf
77 do i3=1,this%nlev
78 read(fileid, iostat=ierr, iomsg=err_message) tempin
79 if (ierr /= 0 ) return
80 this%data(:,i3,i2,i1) = tempin(:)
81 enddo
82 enddo
83 enddo
84 deallocate(tempin)
85 close(fileid)
86
87 end function load
88
89 ! #########################################################################################
90 ! Procedure (type-bound) for setting up interpolation indices between data-grid and
91 ! model-grid.
92 ! #########################################################################################
93 subroutine setup(this, lat, idx1, idx2, idxh)
94 class(ty_h2ophys), intent(in) :: this
95 real(kind_phys), intent(in) :: lat(:)
96 integer, intent(out) :: idx1(:), idx2(:)
97 real(kind_phys), intent(out) :: idxh(:)
98 integer :: i,j
99
100 do j=1,size(lat)
101 idx2(j) = this%nlat + 1
102 do i=1,this%nlat
103 if (lat(j) < this%lat(i)) then
104 idx2(j) = i
105 exit
106 endif
107 enddo
108 idx1(j) = max(idx2(j)-1,1)
109 idx2(j) = min(idx2(j),this%nlat)
110 if (idx2(j) .ne. idx1(j)) then
111 idxh(j) = (lat(j) - this%lat(idx1(j))) / (this%lat(idx2(j)) - this%lat(idx1(j)))
112 else
113 idxh(j) = 1.0
114 endif
115 enddo
116
117 end subroutine setup
118
119 ! #########################################################################################
120 ! Procedure (type-bound) for updating data.
121 ! #########################################################################################
122 subroutine update(this, idx1, idx2, idxh, rjday, idxt1, idxt2, h2opl)
123 class(ty_h2ophys), intent(in) :: this
124 integer, intent(in) :: idx1(:), idx2(:)
125 real(kind_phys), intent(in) :: idxh(:)
126 real(kind_phys), intent(in) :: rjday
127 integer, intent(in) :: idxt1, idxt2
128 real(kind_phys), intent(out) :: h2opl(:,:,:)
129 integer :: nc, l, j, j1, j2
130 real(kind_phys) :: tem, tx1, tx2
131
132 tx1 = (this%time(idxt2) - rjday) / (this%time(idxt2) - this%time(idxt1))
133 tx2 = 1.0 - tx1
134
135 do nc=1,this%ncf
136 do l=1,this%nlev
137 do j=1,size(h2opl(:,1,1))
138 j1 = idx1(j)
139 j2 = idx2(j)
140 tem = 1.0 - idxh(j)
141 h2opl(j,l,nc) = tx1*(tem*this%data(j1,l,nc,idxt1)+idxh(j)*this%data(j2,l,nc,idxt1)) &
142 + tx2*(tem*this%data(j1,l,nc,idxt2)+idxh(j)*this%data(j2,l,nc,idxt2))
143 enddo
144 enddo
145 enddo
146
147 end subroutine update
148
149 ! #########################################################################################
150 ! Procedure (type-bound) for NRL stratospheric h2o photochemistry physics.
151 ! #########################################################################################
152 subroutine run(this, dt, p, h2opltc, h2o, dqv_dt_prd, dqv_dt_qv)
153 class(ty_h2ophys), intent(in) :: this
154 real(kind_phys), intent(in) :: &
155 dt ! Model timestep (sec)
156 real(kind_phys), intent(in), dimension(:,:) :: &
157 p ! Model Pressure (Pa)
158 real(kind_phys), intent(in), dimension(:,:,:) :: &
159 h2opltc ! h2o forcing data
160 real(kind_phys), intent(inout), dimension(:,:) :: &
161 h2o ! h2o concentration (updated)
162 real(kind_phys), intent(inout), dimension(:, :), optional :: &
163 dqv_dt_prd, & ! Net production/loss effect
164 dqv_dt_qv ! water vapor effect
165
166 integer :: nCol, nLev, iCol, iLev, iCf, kmax, kmin, k
167 logical, dimension(size(p,1)) :: flg
168 real(kind_phys) :: pmax, pmin, temp
169 real(kind_phys), dimension(size(p,1)) :: wk1, wk2, wk3, h2oib
170 real(kind_phys), dimension(size(p,1),this%ncf) :: pltc
171 real(kind_phys), parameter :: prsmax=10000.0, pmaxl=log(prsmax)
172
173 ! Dimensions
174 ncol = size(p,1)
175 nlev = size(p,2)
176
177 do ilev=1,nlev
178 pmin = 1.0e10
179 pmax = -1.0e10
180 do icol=1,ncol
181 wk1(icol) = log(p(icol,ilev))
182 pmin = min(wk1(icol), pmin)
183 pmax = max(wk1(icol), pmax)
184 pltc(icol,:) = 0._kind_phys
185 enddo
186 if (pmin < pmaxl) then
187 kmax = 1
188 kmin = 1
189 do k=1,this%nlev-1
190 if (pmin < this%ph2o(k)) kmax = k
191 if (pmax < this%ph2o(k)) kmin = k
192 enddo
193
194 do k=kmin,kmax
195 temp = 1.0 / (this%ph2o(k) - this%ph2o(k+1))
196 do icol=1,ncol
197 flg(icol) = .false.
198 if (wk1(icol) < this%ph2o(k) .and. wk1(icol) >= this%ph2o(k+1)) then
199 flg(icol) = .true.
200 wk2(icol) = (wk1(icol) - this%ph2o(k+1)) * temp
201 wk3(icol) = 1.0 - wk2(icol)
202 endif
203 enddo
204 do icf=1,this%ncf
205 do icol=1,ncol
206 if (flg(icol)) then
207 pltc(icol,icf) = wk2(icol) * h2opltc(icol,k,icf) + wk3(icol) * h2opltc(icol,k+1,icf)
208 endif
209 enddo
210 enddo
211 enddo
212
213 do icf=1,this%ncf
214 do icol=1,ncol
215 if (wk1(icol) < this%ph2o(this%nlev)) then
216 pltc(icol,icf) = h2opltc(icol,this%nlev,icf)
217 endif
218 if (wk1(icol) >= this%ph2o(1)) then
219 pltc(icol,icf) = h2opltc(icol,1,icf)
220 endif
221 enddo
222 enddo
223 endif
224
225 do icol=1,ncol
226 if (p(icol,ilev) < prsmax .and. pltc(icol,2) /= 0.0) then
227 h2oib(icol) = h2o(icol,ilev)
228 temp = 1.0 / pltc(icol,2)
229 h2o(icol,ilev) = (h2oib(icol) + (pltc(icol,1)+pltc(icol,3)*temp)*dt) / (1.0 + temp*dt)
230 endif
231 enddo
232
233 if (present(dqv_dt_prd)) dqv_dt_prd(:, ilev) = pltc(:, 1) * dt
234 if (present(dqv_dt_qv)) then
235 where (pltc(:, 2) .ne. 0)
236 dqv_dt_qv(:, ilev) = (h2o(:, ilev) - pltc(:, 3)) * dt / pltc(:, 2)
237 elsewhere
238 dqv_dt_qv(:, ilev) = 0
239 end where
240 end if
241 enddo
242
243
244 return
245 end subroutine run
246
247end module module_h2ophys
module_h2ophys
Definition module_h2ophys.F90:4
module_h2ophys::ty_h2ophys
Definition module_h2ophys.F90:18