Line data Source code
1 : module wv_saturation
2 :
3 : !--------------------------------------------------------------------!
4 : ! Module Overview: !
5 : ! !
6 : ! This module provides an interface to wv_sat_methods, providing !
7 : ! saturation vapor pressure and related calculations to CAM. !
8 : ! !
9 : ! The original wv_saturation codes were introduced by J. J. Hack, !
10 : ! February 1990. The code has been extensively rewritten since then, !
11 : ! including a total refactoring in Summer 2012. !
12 : ! !
13 : !--------------------------------------------------------------------!
14 : ! Methods: !
15 : ! !
16 : ! Pure water/ice saturation vapor pressures are calculated on the !
17 : ! fly, with the specific method determined by a runtime option. !
18 : ! Mixed phase SVP is interpolated from the internal table, estbl, !
19 : ! which is created during initialization. !
20 : ! !
21 : ! The default method for calculating SVP is determined by a namelist !
22 : ! option, and used whenever svp_water/ice or qsat are called. !
23 : ! !
24 : !--------------------------------------------------------------------!
25 :
26 : use shr_kind_mod, only: r8 => shr_kind_r8
27 : use physconst, only: epsilo, &
28 : latvap, &
29 : latice, &
30 : rh2o, &
31 : cpair, &
32 : tmelt, &
33 : h2otrip
34 :
35 : use wv_sat_methods, only: &
36 : svp_to_qsat => wv_sat_svp_to_qsat, &
37 : svp_to_qsat_vect => wv_sat_svp_to_qsat_vect
38 :
39 : implicit none
40 : private
41 : save
42 :
43 : ! Public interfaces
44 : ! Namelist, initialization, finalization
45 : public wv_sat_readnl
46 : public wv_sat_init
47 : public wv_sat_final
48 :
49 : ! Saturation vapor pressure calculations
50 : public svp_water, svp_water_vect
51 : public svp_ice, svp_ice_vect
52 :
53 : ! Mixed phase (water + ice) saturation vapor pressure table lookup
54 : public estblf
55 :
56 : public svp_to_qsat
57 :
58 : ! Subroutines that return both SVP and humidity
59 : ! Optional arguments do temperature derivatives
60 : interface qsat
61 : module procedure qsat_line
62 : module procedure qsat_vect
63 : module procedure qsat_2D
64 : end interface
65 : public qsat ! Mixed phase
66 : interface qsat_water
67 : module procedure qsat_water_line
68 : module procedure qsat_water_vect
69 : module procedure qsat_water_2D
70 : end interface
71 : public qsat_water ! SVP over water only
72 : interface qsat_ice
73 : module procedure qsat_ice_line
74 : module procedure qsat_ice_vect
75 : module procedure qsat_ice_2D
76 : end interface
77 : public qsat_ice ! SVP over ice only
78 :
79 : ! Wet bulb temperature solver
80 : public :: findsp_vc, findsp
81 :
82 : ! Data
83 :
84 : ! This value is slightly high, but it seems to be the value for the
85 : ! steam point of water originally (and most frequently) used in the
86 : ! Goff & Gratch scheme.
87 : real(r8), parameter :: tboil = 373.16_r8
88 :
89 : ! Table of saturation vapor pressure values (estbl) from tmin to
90 : ! tmax+1 Kelvin, in one degree increments. ttrice defines the
91 : ! transition region, estbl contains a combination of ice & water
92 : ! values.
93 : ! Make these public parameters in case another module wants to see the
94 : ! extent of the table.
95 : real(r8), public, parameter :: tmin = 127.16_r8
96 : real(r8), public, parameter :: tmax = 375.16_r8
97 :
98 : real(r8), parameter :: ttrice = 20.00_r8 ! transition range from es over H2O to es over ice
99 :
100 : integer :: plenest ! length of estbl
101 : real(r8), allocatable :: estbl(:) ! table values of saturation vapor pressure
102 :
103 : real(r8) :: omeps ! 1.0_r8 - epsilo
104 :
105 : real(r8) :: c3 ! parameter used by findsp
106 :
107 : ! Set coefficients for polynomial approximation of difference
108 : ! between saturation vapor press over water and saturation pressure
109 : ! over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial is
110 : ! valid in the range -40 < t < 0 (degrees C).
111 : real(r8) :: pcf(5) = (/ &
112 : 5.04469588506e-01_r8, &
113 : -5.47288442819e+00_r8, &
114 : -3.67471858735e-01_r8, &
115 : -8.95963532403e-03_r8, &
116 : -7.78053686625e-05_r8 /)
117 :
118 : ! --- Degree 6 approximation ---
119 : ! real(r8) :: pcf(6) = (/ &
120 : ! 7.63285250063e-02, &
121 : ! 5.86048427932e+00, &
122 : ! 4.38660831780e-01, &
123 : ! 1.37898276415e-02, &
124 : ! 2.14444472424e-04, &
125 : ! 1.36639103771e-06 /)
126 :
127 : integer, parameter :: VLENS = 128 ! vector length for a GPU kernel
128 :
129 : !$acc declare create (plenest,estbl,omeps,c3,pcf)
130 :
131 : contains
132 :
133 : !---------------------------------------------------------------------
134 : ! ADMINISTRATIVE FUNCTIONS
135 : !---------------------------------------------------------------------
136 :
137 1536 : subroutine wv_sat_readnl(nlfile)
138 : !------------------------------------------------------------------!
139 : ! Purpose: !
140 : ! Get runtime options for wv_saturation. !
141 : !------------------------------------------------------------------!
142 :
143 : use wv_sat_methods, only: wv_sat_get_scheme_idx, &
144 : wv_sat_valid_idx, &
145 : wv_sat_set_default
146 :
147 : use spmd_utils, only: masterproc
148 : use namelist_utils, only: find_group_name
149 : use units, only: getunit, freeunit
150 : use mpishorthand
151 : use cam_abortutils, only: endrun
152 :
153 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
154 :
155 : ! Local variables
156 : integer :: unitn, ierr
157 :
158 : character(len=32) :: wv_sat_scheme = "GoffGratch"
159 :
160 : character(len=*), parameter :: subname = 'wv_sat_readnl'
161 :
162 : namelist /wv_sat_nl/ wv_sat_scheme
163 : !-----------------------------------------------------------------------------
164 :
165 1536 : if (masterproc) then
166 2 : unitn = getunit()
167 2 : open( unitn, file=trim(nlfile), status='old' )
168 2 : call find_group_name(unitn, 'wv_sat_nl', status=ierr)
169 2 : if (ierr == 0) then
170 0 : read(unitn, wv_sat_nl, iostat=ierr)
171 0 : if (ierr /= 0) then
172 0 : call endrun(subname // ':: ERROR reading namelist')
173 0 : return
174 : end if
175 : end if
176 2 : close(unitn)
177 2 : call freeunit(unitn)
178 :
179 : end if
180 :
181 : #ifdef SPMD
182 1536 : call mpibcast(wv_sat_scheme, len(wv_sat_scheme) , mpichar, 0, mpicom)
183 : #endif
184 :
185 1536 : if (.not. wv_sat_set_default(wv_sat_scheme)) then
186 0 : call endrun('wv_sat_readnl :: Invalid wv_sat_scheme.')
187 0 : return
188 : end if
189 :
190 : end subroutine wv_sat_readnl
191 :
192 1536 : subroutine wv_sat_init
193 : !------------------------------------------------------------------!
194 : ! Purpose: !
195 : ! Initialize module (e.g. setting parameters, initializing the !
196 : ! SVP lookup table). !
197 : !------------------------------------------------------------------!
198 :
199 : use wv_sat_methods, only: wv_sat_methods_init, &
200 : wv_sat_get_scheme_idx, &
201 : wv_sat_valid_idx
202 : use spmd_utils, only: masterproc
203 : use cam_logfile, only: iulog
204 : use cam_abortutils, only: endrun
205 : use shr_assert_mod, only: shr_assert_in_domain
206 : use error_messages, only: handle_errmsg
207 :
208 : integer :: status
209 :
210 : ! For wv_sat_methods error reporting.
211 : character(len=256) :: errstring
212 :
213 : ! For generating internal SVP table.
214 : real(r8) :: t ! Temperature
215 : integer :: i ! Increment counter
216 :
217 : ! Precalculated because so frequently used.
218 1536 : omeps = 1.0_r8 - epsilo
219 :
220 : ! Transition range method is only valid for transition temperatures at:
221 : ! -40 deg C < T < 0 deg C
222 : call shr_assert_in_domain(ttrice, ge=0._r8, le=40._r8, varname="ttrice",&
223 1536 : msg="wv_sat_init: Invalid transition temperature range.")
224 :
225 : ! This parameter uses a hardcoded 287.04_r8?
226 1536 : c3 = 287.04_r8*(7.5_r8*log(10._r8))/cpair
227 :
228 : ! Init "methods" module containing actual SVP formulae.
229 :
230 : call wv_sat_methods_init(r8, tmelt, h2otrip, tboil, ttrice, &
231 1536 : epsilo, errstring)
232 :
233 1536 : call handle_errmsg(errstring, subname="wv_sat_methods_init")
234 :
235 : ! Add two to make the table slightly too big, just in case.
236 1536 : plenest = ceiling(tmax-tmin) + 2
237 :
238 : ! Allocate SVP table.
239 1536 : allocate(estbl(plenest), stat=status)
240 1536 : if (status /= 0) then
241 0 : call endrun('wv_sat_init :: ERROR allocating saturation vapor pressure table')
242 0 : return
243 : end if
244 :
245 387072 : do i = 1, plenest
246 387072 : estbl(i) = svp_trans(tmin + real(i-1,r8))
247 : end do
248 :
249 : !$acc update device (plenest,estbl,omeps,c3,pcf)
250 :
251 1536 : if (masterproc) then
252 2 : write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'
253 : end if
254 :
255 : end subroutine wv_sat_init
256 :
257 1536 : subroutine wv_sat_final
258 : !------------------------------------------------------------------!
259 : ! Purpose: !
260 : ! Deallocate global variables in module. !
261 : !------------------------------------------------------------------!
262 : use cam_abortutils, only: endrun
263 :
264 : integer :: status
265 :
266 1536 : if (allocated(estbl)) then
267 :
268 1536 : deallocate(estbl, stat=status)
269 :
270 : if (status /= 0) then
271 : call endrun('wv_sat_final :: ERROR deallocating table')
272 : return
273 : end if
274 :
275 : end if
276 :
277 : end subroutine wv_sat_final
278 :
279 : !---------------------------------------------------------------------
280 : ! DEFAULT SVP FUNCTIONS
281 : !---------------------------------------------------------------------
282 :
283 : ! Compute saturation vapor pressure over water
284 5961927154 : function svp_water(t) result(es)
285 :
286 : use wv_sat_methods, only: &
287 : wv_sat_svp_water
288 :
289 : real(r8), intent(in) :: t ! Temperature (K)
290 : real(r8) :: es ! SVP (Pa)
291 :
292 5961927154 : es = wv_sat_svp_water(t)
293 :
294 5961927154 : end function svp_water
295 :
296 : ! Compute saturation vapor pressure over water
297 750544704 : subroutine svp_water_vect(t, es, vlen)
298 :
299 : use wv_sat_methods, only: &
300 : wv_sat_svp_water_vect
301 :
302 : integer, intent(in) :: vlen
303 : real(r8), intent(in) :: t(vlen) ! Temperature (K)
304 : real(r8), intent(out) :: es(vlen) ! SVP (Pa)
305 :
306 : !$acc data copyin (t) copyout (es)
307 :
308 750544704 : call wv_sat_svp_water_vect(t, es, vlen)
309 :
310 : !$acc end data
311 750544704 : end subroutine svp_water_vect
312 :
313 : ! Compute saturation vapor pressure over ice
314 4691434173 : function svp_ice(t) result(es)
315 :
316 : use wv_sat_methods, only: &
317 : wv_sat_svp_ice
318 :
319 : real(r8), intent(in) :: t ! Temperature (K)
320 : real(r8) :: es ! SVP (Pa)
321 :
322 4691434173 : es = wv_sat_svp_ice(t)
323 :
324 4691434173 : end function svp_ice
325 :
326 : ! Compute saturation vapor pressure over ice
327 750544704 : subroutine svp_ice_vect(t, es, vlen)
328 :
329 : use wv_sat_methods, only: &
330 : wv_sat_svp_ice_vect
331 :
332 : integer, intent(in) :: vlen
333 : real(r8), intent(in) :: t(vlen) ! Temperature (K)
334 : real(r8), intent(out) :: es(vlen) ! SVP (Pa)
335 :
336 : !$acc data copyin(t) copyout(es)
337 :
338 750544704 : call wv_sat_svp_ice_vect(t, es, vlen)
339 :
340 : !$acc end data
341 750544704 : end subroutine svp_ice_vect
342 :
343 : ! Compute saturation vapor pressure with an ice-water transition
344 385536 : function svp_trans(t) result(es)
345 :
346 : use wv_sat_methods, only: &
347 : wv_sat_svp_trans
348 :
349 : real(r8), intent(in) :: t ! Temperature (K)
350 : real(r8) :: es ! SVP (Pa)
351 :
352 385536 : es = wv_sat_svp_trans(t)
353 :
354 385536 : end function svp_trans
355 :
356 : ! Compute saturation vapor pressure with an ice-water transition
357 : subroutine svp_trans_vect(t, es, vlen)
358 :
359 : use wv_sat_methods, only: &
360 : wv_sat_svp_trans_vect
361 :
362 : integer, intent(in) :: vlen
363 : real(r8), intent(in) :: t(vlen) ! Temperature (K)
364 : real(r8), intent(out) :: es(vlen) ! SVP (Pa)
365 :
366 : !$acc data copyin(t) copyout(es)
367 :
368 : call wv_sat_svp_trans_vect(t, es, vlen)
369 :
370 : !$acc end data
371 : end subroutine svp_trans_vect
372 :
373 : !---------------------------------------------------------------------
374 : ! UTILITIES
375 : !---------------------------------------------------------------------
376 :
377 : ! Does linear interpolation from nearest values found
378 : ! in the table (estbl).
379 416923302 : elemental function estblf(t) result(es)
380 :
381 : real(r8), intent(in) :: t ! Temperature
382 : real(r8) :: es ! SVP (Pa)
383 :
384 : integer :: i ! Index for t in the table
385 : real(r8) :: t_tmp ! intermediate temperature for es look-up
386 :
387 : real(r8) :: weight ! Weight for interpolation
388 :
389 416923302 : t_tmp = max(min(t,tmax)-tmin, 0._r8) ! Number of table entries above tmin
390 416923302 : i = int(t_tmp) + 1 ! Corresponding index.
391 416923302 : weight = t_tmp - aint(t_tmp, r8) ! Fractional part of t_tmp (for interpolation).
392 416923302 : es = (1._r8 - weight)*estbl(i) + weight*estbl(i+1)
393 :
394 416923302 : end function estblf
395 :
396 : ! Does linear interpolation from nearest values found
397 : ! in the table (estbl).
398 694819800 : subroutine estblf_vect(t, es, vlen)
399 :
400 : integer, intent(in) :: vlen
401 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
402 : real(r8), dimension(vlen), intent(out) :: es ! SVP (Pa)
403 :
404 : integer :: i ! Index for t in the table
405 : integer :: j
406 : real(r8) :: t_tmp ! intermediate temperature for es look-up
407 :
408 : real(r8) :: weight ! Weight for interpolation
409 :
410 : !$acc data present (t,es)
411 :
412 : !$acc parallel vector_length(VLENS) default(present)
413 : !$acc loop gang vector private(t_tmp,weight,i)
414 11601874800 : do j = 1, vlen
415 10907055000 : t_tmp = max(min(t(j),tmax)-tmin, 0._r8) ! Number of table entries above tmin
416 10907055000 : i = int(t_tmp) + 1 ! Corresponding index.
417 10907055000 : weight = t_tmp - aint(t_tmp, r8) ! Fractional part of t_tmp (for interpolation).
418 11601874800 : es(j) = (1._r8 - weight)*estbl(i) + weight*estbl(i+1)
419 : end do
420 : !$acc end parallel
421 :
422 : !$acc end data
423 694819800 : end subroutine estblf_vect
424 :
425 : ! Get enthalpy based only on temperature
426 : ! and specific humidity.
427 0 : elemental function tq_enthalpy(t, q, hltalt) result(enthalpy)
428 :
429 : real(r8), intent(in) :: t ! Temperature
430 : real(r8), intent(in) :: q ! Specific humidity
431 : real(r8), intent(in) :: hltalt ! Modified hlat for T derivatives
432 :
433 : real(r8) :: enthalpy
434 :
435 0 : enthalpy = cpair * t + hltalt * q
436 :
437 0 : end function tq_enthalpy
438 :
439 : ! Get enthalpy based only on temperature
440 : ! and specific humidity.
441 0 : subroutine tq_enthalpy_vect(t, q, hltalt, enthalpy, vlen)
442 :
443 : integer, intent(in) :: vlen
444 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
445 : real(r8), dimension(vlen), intent(in) :: q ! Specific humidity
446 : real(r8), dimension(vlen), intent(in) :: hltalt ! Modified hlat for T derivatives
447 :
448 : real(r8), dimension(vlen), intent(out) :: enthalpy
449 :
450 : integer :: i
451 :
452 : !$acc data present(t,q,hltalt,enthalpy)
453 :
454 : !$acc parallel vector_length(VLENS) default(present)
455 : !$acc loop gang vector
456 0 : do i = 1, vlen
457 0 : enthalpy(i) = cpair * t(i) + hltalt(i) * q(i)
458 : end do
459 : !$acc end parallel
460 :
461 : !$acc end data
462 0 : end subroutine tq_enthalpy_vect
463 :
464 : !---------------------------------------------------------------------
465 : ! LATENT HEAT OF VAPORIZATION CORRECTIONS
466 : !---------------------------------------------------------------------
467 :
468 0 : elemental subroutine no_ip_hltalt(t, hltalt)
469 : !------------------------------------------------------------------!
470 : ! Purpose: !
471 : ! Calculate latent heat of vaporization of pure liquid water at !
472 : ! a given temperature. !
473 : !------------------------------------------------------------------!
474 :
475 : ! Inputs
476 : real(r8), intent(in) :: t ! Temperature
477 : ! Outputs
478 : real(r8), intent(out) :: hltalt ! Appropriately modified hlat
479 :
480 0 : hltalt = latvap
481 :
482 : ! Account for change of latvap with t above freezing where
483 : ! constant slope is given by -2369 j/(kg c) = cpv - cw
484 0 : if (t >= tmelt) then
485 0 : hltalt = hltalt - 2369.0_r8*(t-tmelt)
486 : end if
487 :
488 0 : end subroutine no_ip_hltalt
489 :
490 375272352 : subroutine no_ip_hltalt_vect(t, hltalt, vlen)
491 : !------------------------------------------------------------------!
492 : ! Purpose: !
493 : ! Calculate latent heat of vaporization of pure liquid water at !
494 : ! a given temperature. !
495 : !------------------------------------------------------------------!
496 :
497 : ! Inputs
498 : integer, intent(in) :: vlen
499 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
500 : ! Outputs
501 : real(r8), dimension(vlen), intent(out) :: hltalt ! Appropriately modified hlat
502 :
503 : integer :: i
504 :
505 : !$acc data present(t,hltalt)
506 :
507 : !$acc parallel vector_length(VLENS) default(present)
508 : !$acc loop gang vector
509 6266175552 : do i = 1, vlen
510 5890903200 : hltalt(i) = latvap
511 : ! Account for change of latvap with t above freezing where
512 : ! constant slope is given by -2369 j/(kg c) = cpv - cw
513 6266175552 : if (t(i) >= tmelt) then
514 1123252109 : hltalt(i) = hltalt(i) - 2369.0_r8*(t(i)-tmelt)
515 : end if
516 : end do
517 : !$acc end parallel
518 :
519 : !$acc end data
520 375272352 : end subroutine no_ip_hltalt_vect
521 :
522 0 : elemental subroutine calc_hltalt(t, hltalt, tterm)
523 : !------------------------------------------------------------------!
524 : ! Purpose: !
525 : ! Calculate latent heat of vaporization of water at a given !
526 : ! temperature, taking into account the ice phase if temperature !
527 : ! is below freezing. !
528 : ! Optional argument also calculates a term used to calculate !
529 : ! d(es)/dT within the water-ice transition range. !
530 : !------------------------------------------------------------------!
531 :
532 : ! Inputs
533 : real(r8), intent(in) :: t ! Temperature
534 : ! Outputs
535 : real(r8), intent(out) :: hltalt ! Appropriately modified hlat
536 : ! Term to account for d(es)/dT in transition region.
537 : real(r8), intent(out), optional :: tterm
538 :
539 : ! Local variables
540 : real(r8) :: tc ! Temperature in degrees C
541 : real(r8) :: weight ! Weight for es transition from water to ice
542 : ! Loop iterator
543 : integer :: i
544 :
545 0 : if (present(tterm)) tterm = 0.0_r8
546 :
547 0 : call no_ip_hltalt(t,hltalt)
548 0 : if (t < tmelt) then
549 : ! Weighting of hlat accounts for transition from water to ice.
550 0 : tc = t - tmelt
551 :
552 0 : if (tc >= -ttrice) then
553 0 : weight = -tc/ttrice
554 :
555 : ! polynomial expression approximates difference between es
556 : ! over water and es over ice from 0 to -ttrice (C) (max of
557 : ! ttrice is 40): required for accurate estimate of es
558 : ! derivative in transition range from ice to water
559 0 : if (present(tterm)) then
560 0 : do i = size(pcf), 1, -1
561 0 : tterm = pcf(i) + tc*tterm
562 : end do
563 0 : tterm = tterm/ttrice
564 : end if
565 :
566 : else
567 : weight = 1.0_r8
568 : end if
569 :
570 0 : hltalt = hltalt + weight*latice
571 :
572 : end if
573 :
574 0 : end subroutine calc_hltalt
575 :
576 0 : subroutine calc_hltalt_vect(t, hltalt, vlen, tterm)
577 : !------------------------------------------------------------------!
578 : ! Purpose: !
579 : ! Calculate latent heat of vaporization of water at a given !
580 : ! temperature, taking into account the ice phase if temperature !
581 : ! is below freezing. !
582 : ! Optional argument also calculates a term used to calculate !
583 : ! d(es)/dT within the water-ice transition range. !
584 : !------------------------------------------------------------------!
585 :
586 : ! Inputs
587 : integer, intent(in) :: vlen
588 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
589 : ! Outputs
590 : real(r8), dimension(vlen), intent(out) :: hltalt ! Appropriately modified hlat
591 : ! Term to account for d(es)/dT in transition region.
592 : real(r8), dimension(vlen), intent(out), optional :: tterm
593 :
594 : ! Local variables
595 : real(r8) :: tc ! Temperature in degrees C
596 : real(r8) :: weight ! Weight for es transition from water to ice
597 : logical :: present_tterm
598 : ! Loop iterator
599 : integer :: i, j, size_pcf
600 :
601 0 : present_tterm = present(tterm)
602 0 : size_pcf = size(pcf)
603 :
604 : !$acc data present(t,hltalt,tterm)
605 :
606 0 : if (present_tterm) then
607 : !$acc parallel vector_length(VLENS) default(present)
608 : !$acc loop gang vector
609 0 : do i = 1, vlen
610 0 : tterm(i) = 0.0_r8
611 : end do
612 : !$acc end parallel
613 : end if
614 :
615 0 : call no_ip_hltalt_vect(t,hltalt,vlen)
616 :
617 : !$acc parallel vector_length(VLENS) default(present)
618 : !$acc loop gang vector private(tc,weight)
619 0 : do j = 1, vlen
620 0 : if (t(j) < tmelt) then
621 : ! Weighting of hlat accounts for transition from water to ice.
622 0 : tc = t(j) - tmelt
623 :
624 0 : if (tc >= -ttrice) then
625 0 : weight = -tc/ttrice
626 :
627 : ! polynomial expression approximates difference between es
628 : ! over water and es over ice from 0 to -ttrice (C) (max of
629 : ! ttrice is 40): required for accurate estimate of es
630 : ! derivative in transition range from ice to water
631 0 : if (present_tterm) then
632 : !$acc loop seq
633 0 : do i = size_pcf, 1, -1
634 0 : tterm(j) = pcf(i) + tc*tterm(j)
635 : end do
636 0 : tterm(j) = tterm(j)/ttrice
637 : end if
638 :
639 : else
640 : weight = 1.0_r8
641 : end if
642 :
643 0 : hltalt(j) = hltalt(j) + weight*latice
644 :
645 : end if
646 : end do
647 : !$acc end parallel
648 :
649 : !$acc end data
650 0 : end subroutine calc_hltalt_vect
651 :
652 : !---------------------------------------------------------------------
653 : ! OPTIONAL OUTPUTS
654 : !---------------------------------------------------------------------
655 :
656 : ! Temperature derivative outputs, for qsat_*
657 0 : subroutine deriv_outputs_line(t, p, es, qs, hltalt, tterm, &
658 : gam, dqsdt)
659 :
660 : ! Inputs
661 : real(r8), intent(in) :: t ! Temperature
662 : real(r8), intent(in) :: p ! Pressure
663 : real(r8), intent(in) :: es ! Saturation vapor pressure
664 : real(r8), intent(in) :: qs ! Saturation specific humidity
665 : real(r8), intent(in) :: hltalt ! Modified latent heat
666 : real(r8), intent(in) :: tterm ! Extra term for d(es)/dT in
667 : ! transition region.
668 :
669 : ! Outputs
670 : real(r8), intent(out), optional :: gam ! (hltalt/cpair)*(d(qs)/dt)
671 : real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt)
672 :
673 : ! Local variables
674 : real(r8) :: desdt ! d(es)/dt
675 : real(r8) :: dqsdt_loc ! local copy of dqsdt
676 :
677 0 : if (qs == 1.0_r8) then
678 : dqsdt_loc = 0._r8
679 : else
680 0 : desdt = hltalt*es/(rh2o*t*t) + tterm
681 0 : dqsdt_loc = qs*p*desdt/(es*(p-omeps*es))
682 : end if
683 :
684 0 : if (present(dqsdt)) dqsdt = dqsdt_loc
685 0 : if (present(gam)) gam = dqsdt_loc * (hltalt/cpair)
686 :
687 0 : end subroutine deriv_outputs_line
688 :
689 : ! Temperature derivative outputs, for qsat_*
690 750544704 : subroutine deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
691 375272352 : gam, dqsdt)
692 :
693 : ! Inputs
694 : integer, intent(in) :: vlen
695 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
696 : real(r8), dimension(vlen), intent(in) :: p ! Pressure
697 : real(r8), dimension(vlen), intent(in) :: es ! Saturation vapor pressure
698 : real(r8), dimension(vlen), intent(in) :: qs ! Saturation specific humidity
699 : real(r8), dimension(vlen), intent(in) :: hltalt ! Modified latent heat
700 : real(r8), dimension(vlen), intent(in) :: tterm ! Extra term for d(es)/dT in
701 : ! transition region.
702 :
703 : ! Outputs
704 : real(r8), dimension(vlen), intent(out), optional :: gam ! (hltalt/cpair)*(d(qs)/dt)
705 : real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt)
706 :
707 : ! Local variables
708 : real(r8) :: desdt ! d(es)/dt
709 : real(r8) :: dqsdt_loc ! local copy of dqsdt
710 : logical :: present_dqsdt, present_gam
711 : integer :: i
712 :
713 375272352 : present_dqsdt = present(dqsdt)
714 375272352 : present_gam = present(gam)
715 :
716 : !$acc data present(t,p,es,qs,hltalt,tterm,gam,dqsdt)
717 :
718 : !$acc parallel vector_length(VLENS) default(present)
719 : !$acc loop gang vector private(dqsdt_loc,desdt)
720 6266175552 : do i = 1, vlen
721 5890903200 : if (qs(i) == 1.0_r8) then
722 : dqsdt_loc = 0._r8
723 : else
724 5782579360 : desdt = hltalt(i)*es(i)/(rh2o*t(i)*t(i)) + tterm(i)
725 5782579360 : dqsdt_loc = qs(i)*p(i)*desdt/(es(i)*(p(i)-omeps*es(i)))
726 : end if
727 :
728 5890903200 : if (present_dqsdt) dqsdt(i) = dqsdt_loc
729 6266175552 : if (present_gam) gam(i) = dqsdt_loc * (hltalt(i)/cpair)
730 : end do
731 : !$acc end parallel
732 :
733 : !$acc end data
734 375272352 : end subroutine deriv_outputs_vect
735 :
736 : !---------------------------------------------------------------------
737 : ! QSAT (SPECIFIC HUMIDITY) PROCEDURES
738 : !---------------------------------------------------------------------
739 :
740 416923302 : subroutine qsat_line(t, p, es, qs, gam, dqsdt, enthalpy)
741 : !------------------------------------------------------------------!
742 : ! Purpose: !
743 : ! Look up and return saturation vapor pressure from precomputed !
744 : ! table, then calculate and return saturation specific humidity. !
745 : ! Optionally return various temperature derivatives or enthalpy !
746 : ! at saturation. !
747 : !------------------------------------------------------------------!
748 :
749 : ! Inputs
750 : real(r8), intent(in) :: t ! Temperature
751 : real(r8), intent(in) :: p ! Pressure
752 : ! Outputs
753 : real(r8), intent(out) :: es ! Saturation vapor pressure
754 : real(r8), intent(out) :: qs ! Saturation specific humidity
755 :
756 : real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
757 : real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt)
758 : real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
759 :
760 : ! Local variables
761 : real(r8) :: hltalt ! Modified latent heat for T derivatives
762 : real(r8) :: tterm ! Account for d(es)/dT in transition region
763 :
764 416923302 : es = estblf(t)
765 :
766 416923302 : qs = svp_to_qsat(es, p)
767 :
768 : ! Ensures returned es is consistent with limiters on qs.
769 416923302 : es = min(es, p)
770 :
771 : ! Calculate optional arguments.
772 416923302 : if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then
773 :
774 : ! "generalized" analytic expression for t derivative of es
775 : ! accurate to within 1 percent for 173.16 < t < 373.16
776 0 : call calc_hltalt(t, hltalt, tterm)
777 :
778 0 : if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt)
779 :
780 : call deriv_outputs_line(t, p, es, qs, hltalt, tterm, &
781 0 : gam=gam, dqsdt=dqsdt)
782 :
783 : end if
784 :
785 416923302 : end subroutine qsat_line
786 :
787 694819800 : subroutine qsat_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy)
788 : !------------------------------------------------------------------!
789 : ! Purpose: !
790 : ! Look up and return saturation vapor pressure from precomputed !
791 : ! table, then calculate and return saturation specific humidity. !
792 : ! Optionally return various temperature derivatives or enthalpy !
793 : ! at saturation. !
794 : !------------------------------------------------------------------!
795 :
796 : ! Inputs
797 : integer, intent(in) :: vlen
798 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
799 : real(r8), dimension(vlen), intent(in) :: p ! Pressure
800 : ! Outputs
801 : real(r8), dimension(vlen), intent(out) :: es ! Saturation vapor pressure
802 : real(r8), dimension(vlen), intent(out) :: qs ! Saturation specific humidity
803 :
804 : real(r8), dimension(vlen), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
805 : real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt)
806 : real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
807 :
808 : ! Local variables
809 1389639600 : real(r8), dimension(vlen) :: hltalt ! Modified latent heat for T derivatives
810 1389639600 : real(r8), dimension(vlen) :: tterm ! Account for d(es)/dT in transition region
811 : integer :: i
812 : logical :: present_gam, present_dqsdt, present_enthalpy
813 :
814 694819800 : present_gam = present(gam)
815 694819800 : present_dqsdt = present(dqsdt)
816 694819800 : present_enthalpy = present(enthalpy)
817 :
818 : !$acc data copyin (t,p) &
819 : !$acc copyout (es,qs,gam,dqsdt,enthalpy) &
820 : !$acc create (hltalt,tterm)
821 :
822 694819800 : call estblf_vect(t, es, vlen)
823 :
824 694819800 : call svp_to_qsat_vect(es, p, qs, vlen)
825 :
826 : ! Ensures returned es is consistent with limiters on qs.
827 :
828 : !$acc parallel vector_length(VLENS) default(present)
829 : !$acc loop gang vector
830 11601874800 : do i = 1, vlen
831 11601874800 : es(i) = min(es(i), p(i))
832 : end do
833 : !$acc end parallel
834 :
835 : ! Calculate optional arguments.
836 694819800 : if (present_gam .or. present_dqsdt .or. present_enthalpy) then
837 :
838 : ! "generalized" analytic expression for t derivative of es
839 : ! accurate to within 1 percent for 173.16 < t < 373.16
840 0 : call calc_hltalt_vect(t, hltalt, vlen, tterm)
841 :
842 0 : if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
843 :
844 : call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
845 0 : gam=gam, dqsdt=dqsdt)
846 :
847 : end if
848 :
849 : !$acc end data
850 694819800 : end subroutine qsat_vect
851 :
852 0 : subroutine qsat_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy)
853 : !------------------------------------------------------------------!
854 : ! Purpose: !
855 : ! Look up and return saturation vapor pressure from precomputed !
856 : ! table, then calculate and return saturation specific humidity. !
857 : ! Optionally return various temperature derivatives or enthalpy !
858 : ! at saturation. !
859 : !------------------------------------------------------------------!
860 :
861 : ! Inputs
862 : integer, intent(in) :: dim1, dim2
863 : real(r8), dimension(dim1,dim2), intent(in) :: t ! Temperature
864 : real(r8), dimension(dim1,dim2), intent(in) :: p ! Pressure
865 : ! Outputs
866 : real(r8), dimension(dim1,dim2), intent(out) :: es ! Saturation vapor pressure
867 : real(r8), dimension(dim1,dim2), intent(out) :: qs ! Saturation specific humidity
868 :
869 : real(r8), dimension(dim1,dim2), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
870 : real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt ! (d(qs)/dt)
871 : real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
872 :
873 : ! Local variables
874 0 : real(r8), dimension(dim1,dim2) :: hltalt ! Modified latent heat for T derivatives
875 0 : real(r8), dimension(dim1,dim2) :: tterm ! Account for d(es)/dT in transition region
876 : integer :: i, k, vlen
877 : logical :: present_gam, present_dqsdt, present_enthalpy
878 :
879 0 : vlen = dim1 * dim2
880 0 : present_gam = present(gam)
881 0 : present_dqsdt = present(dqsdt)
882 0 : present_enthalpy = present(enthalpy)
883 :
884 : !$acc data copyin (t,p) &
885 : !$acc copyout (es,qs,gam,dqsdt,enthalpy) &
886 : !$acc create (hltalt,tterm)
887 :
888 0 : call estblf_vect(t, es, vlen)
889 :
890 0 : call svp_to_qsat_vect(es, p, qs, vlen)
891 :
892 : ! Ensures returned es is consistent with limiters on qs.
893 :
894 : !$acc parallel vector_length(VLENS) default(present)
895 : !$acc loop gang vector collapse(2)
896 0 : do k = 1, dim2
897 0 : do i = 1, dim1
898 0 : es(i,k) = min(es(i,k), p(i,k))
899 : end do
900 : end do
901 : !$acc end parallel
902 :
903 : ! Calculate optional arguments.
904 0 : if (present_gam .or. present_dqsdt .or. present_enthalpy) then
905 :
906 : ! "generalized" analytic expression for t derivative of es
907 : ! accurate to within 1 percent for 173.16 < t < 373.16
908 0 : call calc_hltalt_vect(t, hltalt, vlen, tterm)
909 :
910 0 : if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
911 :
912 : call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
913 0 : gam=gam, dqsdt=dqsdt)
914 :
915 : end if
916 :
917 : !$acc end data
918 0 : end subroutine qsat_2D
919 :
920 29058269913 : subroutine qsat_water_line(t, p, es, qs, gam, dqsdt, enthalpy)
921 : !------------------------------------------------------------------!
922 : ! Purpose: !
923 : ! Calculate SVP over water at a given temperature, and then !
924 : ! calculate and return saturation specific humidity. !
925 : ! Optionally return various temperature derivatives or enthalpy !
926 : ! at saturation. !
927 : !------------------------------------------------------------------!
928 :
929 : use wv_sat_methods, only: wv_sat_qsat_water
930 :
931 : ! Inputs
932 : real(r8), intent(in) :: t ! Temperature
933 : real(r8), intent(in) :: p ! Pressure
934 : ! Outputs
935 : real(r8), intent(out) :: es ! Saturation vapor pressure
936 : real(r8), intent(out) :: qs ! Saturation specific humidity
937 :
938 : real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
939 : real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt)
940 : real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
941 :
942 : ! Local variables
943 : real(r8) :: hltalt ! Modified latent heat for T derivatives
944 :
945 29058269913 : call wv_sat_qsat_water(t, p, es, qs)
946 :
947 29058269913 : if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then
948 :
949 : ! "generalized" analytic expression for t derivative of es
950 : ! accurate to within 1 percent for 173.16 < t < 373.16
951 0 : call no_ip_hltalt(t, hltalt)
952 :
953 0 : if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt)
954 :
955 : ! For pure water/ice transition term is 0.
956 : call deriv_outputs_line(t, p, es, qs, hltalt, 0._r8, &
957 0 : gam=gam, dqsdt=dqsdt)
958 :
959 : end if
960 :
961 29058269913 : end subroutine qsat_water_line
962 :
963 889326000 : subroutine qsat_water_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy)
964 : !------------------------------------------------------------------!
965 : ! Purpose: !
966 : ! Calculate SVP over water at a given temperature, and then !
967 : ! calculate and return saturation specific humidity. !
968 : ! Optionally return various temperature derivatives or enthalpy !
969 : ! at saturation. !
970 : !------------------------------------------------------------------!
971 :
972 : use wv_sat_methods, only: wv_sat_qsat_water_vect
973 :
974 : ! Inputs
975 : integer, intent(in) :: vlen
976 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
977 : real(r8), dimension(vlen), intent(in) :: p ! Pressure
978 : ! Outputs
979 : real(r8), dimension(vlen), intent(out) :: es ! Saturation vapor pressure
980 : real(r8), dimension(vlen), intent(out) :: qs ! Saturation specific humidity
981 :
982 : real(r8), dimension(vlen), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
983 : real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt)
984 : real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
985 :
986 : ! Local variables
987 1778652000 : real(r8), dimension(vlen) :: hltalt ! Modified latent heat for T derivatives
988 1778652000 : real(r8), dimension(vlen) :: tterm
989 : integer :: i
990 : logical :: present_gam, present_dqsdt, present_enthalpy
991 :
992 889326000 : present_gam = present(gam)
993 889326000 : present_dqsdt = present(dqsdt)
994 889326000 : present_enthalpy = present(enthalpy)
995 :
996 : !$acc data copyin (t,p) &
997 : !$acc copyout (es,qs,gam,dqsdt,enthalpy) &
998 : !$acc create (tterm,hltalt)
999 :
1000 : !$acc parallel vector_length(VLENS) default(present)
1001 : !$acc loop gang vector
1002 14849676000 : do i = 1, vlen
1003 14849676000 : tterm(i) = 0._r8
1004 : end do
1005 : !$acc end parallel
1006 :
1007 889326000 : call wv_sat_qsat_water_vect(t, p, es, qs, vlen)
1008 :
1009 889326000 : if (present_gam .or. present_dqsdt .or. present_enthalpy) then
1010 :
1011 : ! "generalized" analytic expression for t derivative of es
1012 : ! accurate to within 1 percent for 173.16 < t < 373.16
1013 375272352 : call no_ip_hltalt_vect(t, hltalt, vlen)
1014 :
1015 375272352 : if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
1016 :
1017 : ! For pure water/ice transition term is 0.
1018 : call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
1019 375272352 : gam=gam, dqsdt=dqsdt)
1020 :
1021 : end if
1022 :
1023 : !$acc end data
1024 889326000 : end subroutine qsat_water_vect
1025 :
1026 0 : subroutine qsat_water_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy)
1027 : !------------------------------------------------------------------!
1028 : ! Purpose: !
1029 : ! Calculate SVP over water at a given temperature, and then !
1030 : ! calculate and return saturation specific humidity. !
1031 : ! Optionally return various temperature derivatives or enthalpy !
1032 : ! at saturation. !
1033 : !------------------------------------------------------------------!
1034 :
1035 : use wv_sat_methods, only: wv_sat_qsat_water_vect
1036 :
1037 : ! Inputs
1038 : integer, intent(in) :: dim1, dim2
1039 : real(r8), dimension(dim1,dim2), intent(in) :: t ! Temperature
1040 : real(r8), dimension(dim1,dim2), intent(in) :: p ! Pressure
1041 : ! Outputs
1042 : real(r8), dimension(dim1,dim2), intent(out) :: es ! Saturation vapor pressure
1043 : real(r8), dimension(dim1,dim2), intent(out) :: qs ! Saturation specific humidity
1044 :
1045 : real(r8), dimension(dim1,dim2), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
1046 : real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt ! (d(qs)/dt)
1047 : real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
1048 :
1049 : ! Local variables
1050 0 : real(r8), dimension(dim1,dim2) :: hltalt ! Modified latent heat for T derivatives
1051 0 : real(r8), dimension(dim1,dim2) :: tterm
1052 : integer :: i, k, vlen
1053 : logical :: present_gam, present_dqsdt, present_enthalpy
1054 :
1055 0 : vlen = dim1 * dim2
1056 0 : present_gam = present(gam)
1057 0 : present_dqsdt = present(dqsdt)
1058 0 : present_enthalpy = present(enthalpy)
1059 :
1060 : !$acc data copyin (t,p) &
1061 : !$acc copyout (es,qs,gam,dqsdt,enthalpy) &
1062 : !$acc create (hltalt,tterm)
1063 :
1064 : !$acc parallel vector_length(VLENS) default(present)
1065 : !$acc loop gang vector collapse(2)
1066 0 : do k = 1, dim2
1067 0 : do i = 1, dim1
1068 0 : tterm(i,k) = 0._r8
1069 : end do
1070 : end do
1071 : !$acc end parallel
1072 :
1073 0 : call wv_sat_qsat_water_vect(t, p, es, qs, vlen)
1074 :
1075 0 : if (present_gam .or. present_dqsdt .or. present_enthalpy) then
1076 :
1077 : ! "generalized" analytic expression for t derivative of es
1078 : ! accurate to within 1 percent for 173.16 < t < 373.16
1079 0 : call no_ip_hltalt_vect(t, hltalt, vlen)
1080 :
1081 0 : if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
1082 :
1083 : ! For pure water/ice transition term is 0.
1084 : call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
1085 0 : gam=gam, dqsdt=dqsdt)
1086 :
1087 : end if
1088 :
1089 : !$acc end data
1090 0 : end subroutine qsat_water_2D
1091 :
1092 0 : subroutine qsat_ice_line(t, p, es, qs, gam, dqsdt, enthalpy)
1093 : !------------------------------------------------------------------!
1094 : ! Purpose: !
1095 : ! Calculate SVP over ice at a given temperature, and then !
1096 : ! calculate and return saturation specific humidity. !
1097 : ! Optionally return various temperature derivatives or enthalpy !
1098 : ! at saturation. !
1099 : !------------------------------------------------------------------!
1100 :
1101 : use wv_sat_methods, only: wv_sat_qsat_ice
1102 :
1103 : ! Inputs
1104 : real(r8), intent(in) :: t ! Temperature
1105 : real(r8), intent(in) :: p ! Pressure
1106 : ! Outputs
1107 : real(r8), intent(out) :: es ! Saturation vapor pressure
1108 : real(r8), intent(out) :: qs ! Saturation specific humidity
1109 :
1110 : real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
1111 : real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt)
1112 : real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
1113 :
1114 : ! Local variables
1115 : real(r8) :: hltalt ! Modified latent heat for T derivatives
1116 :
1117 0 : call wv_sat_qsat_ice(t, p, es, qs)
1118 :
1119 0 : if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then
1120 :
1121 : ! For pure ice, just add latent heats.
1122 0 : hltalt = latvap + latice
1123 :
1124 0 : if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt)
1125 :
1126 : ! For pure water/ice transition term is 0.
1127 : call deriv_outputs_line(t, p, es, qs, hltalt, 0._r8, &
1128 0 : gam=gam, dqsdt=dqsdt)
1129 :
1130 : end if
1131 :
1132 0 : end subroutine qsat_ice_line
1133 :
1134 0 : subroutine qsat_ice_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy)
1135 : !------------------------------------------------------------------!
1136 : ! Purpose: !
1137 : ! Calculate SVP over ice at a given temperature, and then !
1138 : ! calculate and return saturation specific humidity. !
1139 : ! Optionally return various temperature derivatives or enthalpy !
1140 : ! at saturation. !
1141 : !------------------------------------------------------------------!
1142 :
1143 : use wv_sat_methods, only: wv_sat_qsat_ice_vect
1144 :
1145 : ! Inputs
1146 : integer, intent(in) :: vlen
1147 : real(r8), dimension(vlen), intent(in) :: t ! Temperature
1148 : real(r8), dimension(vlen), intent(in) :: p ! Pressure
1149 : ! Outputs
1150 : real(r8), dimension(vlen), intent(out) :: es ! Saturation vapor pressure
1151 : real(r8), dimension(vlen), intent(out) :: qs ! Saturation specific humidity
1152 :
1153 : real(r8), dimension(vlen), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
1154 : real(r8), dimension(vlen), intent(out), optional :: dqsdt ! (d(qs)/dt)
1155 : real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
1156 :
1157 : ! Local variables
1158 0 : real(r8), dimension(vlen) :: hltalt ! Modified latent heat for T derivatives
1159 0 : real(r8), dimension(vlen) :: tterm
1160 : integer :: i
1161 : logical :: present_gam, present_dqsdt, present_enthalpy
1162 :
1163 0 : present_gam = present(gam)
1164 0 : present_dqsdt = present(dqsdt)
1165 0 : present_enthalpy = present(enthalpy)
1166 :
1167 : !$acc data copyin (t,p) &
1168 : !$acc copyout (es,qs,gam,dqsdt,enthalpy) &
1169 : !$acc create (hltalt,tterm)
1170 :
1171 : !$acc parallel vector_length(VLENS) default(present)
1172 : !$acc loop gang vector
1173 0 : do i = 1, vlen
1174 0 : tterm(i) = 0._r8
1175 : end do
1176 : !$acc end parallel
1177 :
1178 0 : call wv_sat_qsat_ice_vect(t, p, es, qs, vlen)
1179 :
1180 0 : if (present_gam .or. present_dqsdt .or. present_enthalpy) then
1181 :
1182 : !$acc parallel vector_length(VLENS) default(present)
1183 : !$acc loop gang vector
1184 0 : do i = 1, vlen
1185 : ! For pure ice, just add latent heats.
1186 0 : hltalt(i) = latvap + latice
1187 : end do
1188 : !$acc end parallel
1189 :
1190 0 : if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
1191 :
1192 : ! For pure water/ice transition term is 0.
1193 : call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
1194 0 : gam=gam, dqsdt=dqsdt)
1195 :
1196 : end if
1197 :
1198 : !$acc end data
1199 0 : end subroutine qsat_ice_vect
1200 :
1201 0 : subroutine qsat_ice_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy)
1202 : !------------------------------------------------------------------!
1203 : ! Purpose: !
1204 : ! Calculate SVP over ice at a given temperature, and then !
1205 : ! calculate and return saturation specific humidity. !
1206 : ! Optionally return various temperature derivatives or enthalpy !
1207 : ! at saturation. !
1208 : !------------------------------------------------------------------!
1209 :
1210 : use wv_sat_methods, only: wv_sat_qsat_ice_vect
1211 :
1212 : ! Inputs
1213 : integer, intent(in) :: dim1, dim2
1214 : real(r8), dimension(dim1,dim2), intent(in) :: t ! Temperature
1215 : real(r8), dimension(dim1,dim2), intent(in) :: p ! Pressure
1216 : ! Outputs
1217 : real(r8), dimension(dim1,dim2), intent(out) :: es ! Saturation vapor pressure
1218 : real(r8), dimension(dim1,dim2), intent(out) :: qs ! Saturation specific humidity
1219 :
1220 : real(r8), dimension(dim1,dim2), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt)
1221 : real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt ! (d(qs)/dt)
1222 : real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
1223 :
1224 : ! Local variables
1225 0 : real(r8), dimension(dim1,dim2) :: hltalt ! Modified latent heat for T derivatives
1226 0 : real(r8), dimension(dim1,dim2) :: tterm
1227 : integer :: i, k, vlen
1228 : logical :: present_gam, present_dqsdt, present_enthalpy
1229 :
1230 0 : vlen = dim1 * dim2
1231 0 : present_gam = present(gam)
1232 0 : present_dqsdt = present(dqsdt)
1233 0 : present_enthalpy = present(enthalpy)
1234 :
1235 : !$acc data copyin (t,p) &
1236 : !$acc copyout (es,qs,gam,dqsdt,enthalpy) &
1237 : !$acc create (hltalt,tterm)
1238 :
1239 : !$acc parallel vector_length(VLENS) default(present)
1240 : !$acc loop gang vector collapse(2)
1241 0 : do k = 1, dim2
1242 0 : do i = 1, dim1
1243 0 : tterm(i,k) = 0._r8
1244 : end do
1245 : end do
1246 : !$acc end parallel
1247 :
1248 0 : call wv_sat_qsat_ice_vect(t, p, es, qs, vlen)
1249 :
1250 0 : if (present_gam .or. present_dqsdt .or. present_enthalpy) then
1251 :
1252 : !$acc parallel vector_length(VLENS) default(present)
1253 : !$acc loop gang vector collapse(2)
1254 0 : do k = 1, dim2
1255 0 : do i = 1, dim1
1256 : ! For pure ice, just add latent heats.
1257 0 : hltalt(i,k) = latvap + latice
1258 : end do
1259 : end do
1260 : !$acc end parallel
1261 :
1262 0 : if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
1263 :
1264 : ! For pure water/ice transition term is 0.
1265 : call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
1266 0 : gam=gam, dqsdt=dqsdt)
1267 :
1268 : end if
1269 :
1270 : !$acc end data
1271 0 : end subroutine qsat_ice_2D
1272 :
1273 : !---------------------------------------------------------------------
1274 : ! FINDSP (WET BULB TEMPERATURE) PROCEDURES
1275 : !---------------------------------------------------------------------
1276 :
1277 0 : subroutine findsp_vc(q, t, p, use_ice, tsp, qsp)
1278 :
1279 : use cam_logfile, only: iulog
1280 : use cam_abortutils, only: endrun
1281 :
1282 : ! Wrapper for findsp which is 1D and handles the output status.
1283 : ! Changing findsp to elemental restricted debugging output.
1284 : ! If that output is needed again, it's preferable *not* to copy findsp,
1285 : ! but to change the existing version.
1286 :
1287 : ! input arguments
1288 : real(r8), intent(in) :: q(:) ! water vapor (kg/kg)
1289 : real(r8), intent(in) :: t(:) ! temperature (K)
1290 : real(r8), intent(in) :: p(:) ! pressure (Pa)
1291 : logical, intent(in) :: use_ice ! flag to include ice phase in calculations
1292 :
1293 : ! output arguments
1294 : real(r8), intent(out) :: tsp(:) ! saturation temp (K)
1295 : real(r8), intent(out) :: qsp(:) ! saturation mixing ratio (kg/kg)
1296 :
1297 0 : integer :: status(size(q)) ! flag representing state of output
1298 : ! 0 => Successful convergence
1299 : ! 1 => No calculation done: pressure or specific
1300 : ! humidity not within usable range
1301 : ! 2 => Run failed to converge
1302 : ! 4 => Temperature fell below minimum
1303 : ! 8 => Enthalpy not conserved
1304 :
1305 : integer :: n, i
1306 :
1307 0 : n = size(q)
1308 :
1309 : ! Currently, only 2 and 8 seem to be treated as fatal errors.
1310 0 : do i = 1,n
1311 0 : call findsp(q(i), t(i), p(i), use_ice, tsp(i), qsp(i), status(i))
1312 0 : if (status(i) == 2) then
1313 0 : write(iulog,*) ' findsp not converging at i = ', i
1314 0 : write(iulog,*) ' t, q, p ', t(i), q(i), p(i)
1315 0 : write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i)
1316 0 : call endrun ('wv_saturation::FINDSP -- not converging')
1317 0 : else if (status(i) == 8) then
1318 0 : write(iulog,*) ' the enthalpy is not conserved at i = ', i
1319 0 : write(iulog,*) ' t, q, p ', t(i), q(i), p(i)
1320 0 : write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i)
1321 0 : call endrun ('wv_saturation::FINDSP -- enthalpy is not conserved')
1322 : endif
1323 : end do
1324 :
1325 0 : end subroutine findsp_vc
1326 :
1327 0 : subroutine findsp (q, t, p, use_ice, tsp, qsp, status)
1328 : !-----------------------------------------------------------------------
1329 : !
1330 : ! Purpose:
1331 : ! find the wet bulb temperature for a given t and q
1332 : ! in a longitude height section
1333 : ! wet bulb temp is the temperature and spec humidity that is
1334 : ! just saturated and has the same enthalpy
1335 : ! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
1336 : ! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
1337 : !
1338 : ! Method:
1339 : ! a Newton method is used
1340 : ! first guess uses an algorithm provided by John Petch from the UKMO
1341 : ! we exclude points where the physical situation is unrealistic
1342 : ! e.g. where the temperature is outside the range of validity for the
1343 : ! saturation vapor pressure, or where the water vapor pressure
1344 : ! exceeds the ambient pressure, or the saturation specific humidity is
1345 : ! unrealistic
1346 : !
1347 : ! Author: P. Rasch
1348 : !
1349 : !-----------------------------------------------------------------------
1350 : !
1351 : ! input arguments
1352 : !
1353 :
1354 : real(r8), intent(in) :: q ! water vapor (kg/kg)
1355 : real(r8), intent(in) :: t ! temperature (K)
1356 : real(r8), intent(in) :: p ! pressure (Pa)
1357 : logical, intent(in) :: use_ice ! flag to include ice phase in calculations
1358 : !
1359 : ! output arguments
1360 : !
1361 : real(r8), intent(out) :: tsp ! saturation temp (K)
1362 : real(r8), intent(out) :: qsp ! saturation mixing ratio (kg/kg)
1363 : integer, intent(out) :: status ! flag representing state of output
1364 : ! 0 => Successful convergence
1365 : ! 1 => No calculation done: pressure or specific
1366 : ! humidity not within usable range
1367 : ! 2 => Run failed to converge
1368 : ! 4 => Temperature fell below minimum
1369 : ! 8 => Enthalpy not conserved
1370 : !
1371 : ! local variables
1372 : !
1373 : integer, parameter :: iter = 8 ! max number of times to iterate the calculation
1374 : integer :: l ! iterator
1375 :
1376 : real(r8) es ! sat. vapor pressure
1377 : real(r8) gam ! change in sat spec. hum. wrt temperature (times hltalt/cpair)
1378 : real(r8) dgdt ! work variable
1379 : real(r8) g ! work variable
1380 : real(r8) hltalt ! lat. heat. of vap.
1381 : real(r8) qs ! spec. hum. of water vapor
1382 :
1383 : ! work variables
1384 : real(r8) t1, q1, dt, dq
1385 : real(r8) qvd
1386 : real(r8) r1b, c1, c2
1387 : real(r8), parameter :: dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
1388 : real(r8), parameter :: dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
1389 : real(r8) enin, enout
1390 :
1391 : ! Saturation specific humidity at this temperature
1392 0 : if (use_ice) then
1393 0 : call qsat(t, p, es, qs)
1394 : else
1395 0 : call qsat_water(t, p, es, qs)
1396 : end if
1397 :
1398 : ! make sure a meaningful calculation is possible
1399 : if (p <= 5._r8*es .or. qs <= 0._r8 .or. qs >= 0.5_r8 &
1400 0 : .or. t < tmin .or. t > tmax) then
1401 0 : status = 1
1402 : ! Keep initial parameters when conditions aren't suitable
1403 0 : tsp = t
1404 0 : qsp = q
1405 0 : enin = 1._r8
1406 0 : enout = 1._r8
1407 :
1408 0 : return
1409 : end if
1410 :
1411 : ! Prepare to iterate
1412 0 : status = 2
1413 :
1414 : ! Get initial enthalpy
1415 0 : if (use_ice) then
1416 0 : call calc_hltalt(t,hltalt)
1417 : else
1418 0 : call no_ip_hltalt(t,hltalt)
1419 : end if
1420 0 : enin = tq_enthalpy(t, q, hltalt)
1421 :
1422 : ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
1423 0 : c1 = hltalt*c3
1424 0 : c2 = (t + 36._r8)**2
1425 0 : r1b = c2/(c2 + c1*qs)
1426 0 : qvd = r1b * (q - qs)
1427 0 : tsp = t + ((hltalt/cpair)*qvd)
1428 :
1429 : ! Generate qsp, gam, and enout from tsp.
1430 0 : if (use_ice) then
1431 0 : call qsat(tsp, p, es, qsp, gam=gam, enthalpy=enout)
1432 : else
1433 0 : call qsat_water(tsp, p, es, qsp, gam=gam, enthalpy=enout)
1434 : end if
1435 :
1436 : ! iterate on first guess
1437 0 : do l = 1, iter
1438 :
1439 0 : g = enin - enout
1440 0 : dgdt = -cpair * (1 + gam)
1441 :
1442 : ! New tsp
1443 0 : t1 = tsp - g/dgdt
1444 0 : dt = abs(t1 - tsp)/t1
1445 0 : tsp = t1
1446 :
1447 : ! bail out if past end of temperature range
1448 0 : if ( tsp < tmin ) then
1449 0 : tsp = tmin
1450 : ! Get latent heat and set qsp to a value
1451 : ! that preserves enthalpy.
1452 0 : if (use_ice) then
1453 0 : call calc_hltalt(tsp,hltalt)
1454 : else
1455 0 : call no_ip_hltalt(tsp,hltalt)
1456 : end if
1457 0 : qsp = (enin - cpair*tsp)/hltalt
1458 0 : enout = tq_enthalpy(tsp, qsp, hltalt)
1459 0 : status = 4
1460 0 : exit
1461 : end if
1462 :
1463 : ! Re-generate qsp, gam, and enout from new tsp.
1464 0 : if (use_ice) then
1465 0 : call qsat(tsp, p, es, q1, gam=gam, enthalpy=enout)
1466 : else
1467 0 : call qsat_water(tsp, p, es, q1, gam=gam, enthalpy=enout)
1468 : end if
1469 0 : dq = abs(q1 - qsp)/max(q1,1.e-12_r8)
1470 0 : qsp = q1
1471 :
1472 : ! if converged at this point, exclude it from more iterations
1473 0 : if (dt < dttol .and. dq < dqtol) then
1474 0 : status = 0
1475 0 : exit
1476 : endif
1477 : end do
1478 :
1479 : ! Test for enthalpy conservation
1480 0 : if (abs((enin-enout)/(enin+enout)) > 1.e-4_r8) status = 8
1481 :
1482 : end subroutine findsp
1483 :
1484 : end module wv_saturation
|
