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