CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
wv_saturation.F
1
4!
5! $Id: wv_saturation.F90,v 1.1.12.2.10.1 2014-04-14 16:04:56 dbarahon Exp $
6!
11#ifdef GEOS5
12 use mapl_constantsmod, kp => mapl_r8
13#endif
14#ifdef NEMS_GSM
15 use funcphys, only : fpvsl, fpvsi, fpvs ! saturation vapor pressure for water & ice
16 use machine, only : kp => kind_phys
17#endif
18
19!++jtb (comm out)
20
21
22
23!--jtb
24
25 implicit none
26 private
27 save
28!
29! Public interfaces
30!
31 public gestbl
32 public estblf
33 public aqsat
34 public aqsatd
35 public vqsatd
36 public fqsatd
37 public qsat_water
38 public vqsat_water
39 public qsat_ice
40 public vqsat_ice
41 public vqsatd_water
42 public aqsat_water
43 public vqsatd2_water
46 public vqsatd2
47 public vqsatd2_single
48 public polysvp
49!
50! Data used by cldwat
51!
52 public hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice
53!
54! Data
55!
56 integer plenest
57 parameter(plenest=250)
58!
59! Table of saturation vapor pressure values es from tmin degrees
60! to tmax+1 degrees k in one degree increments. ttrice defines the
61! transition region where es is a combination of ice & water values
62!
63 real(kp) estbl(plenest)
64 real(kp) tmin
65 real(kp) tmax
66 real(kp) ttrice
67 real(kp) pcf(6)
68 real(kp) epsqs
69 real(kp) rgasv
70 real(kp) hlatf
71 real(kp) hlatv
72 real(kp) cp
73 real(kp) tmelt
74 logical icephs
75
76 integer, parameter :: iulog=6
77
78 contains
79
80 real(kp) function estblf( td )
81!
82! Saturation vapor pressure table lookup
83!
84 real(kp), intent(in) :: td
85!
86 real(kp) :: e
87 real(kp) :: ai
88 integer :: i
89!
90 e = max(min(td,tmax),tmin)
91 i = int(e-tmin)+1
92 ai = aint(e-tmin)
93 estblf = (tmin+ai-e+1._kp)* estbl(i)-(tmin+ai-e)* estbl(i+1)
94 end function estblf
95
106 subroutine gestbl(tmn ,tmx ,trice ,ip ,epsil , latvap ,latice , &
107 & rh2o ,cpair ,tmeltx )
108!------------------------------Arguments--------------------------------
109!
110! Input arguments
111!
112 real(kp), intent(in) :: tmn
113 real(kp), intent(in) :: tmx
114 real(kp), intent(in) :: epsil
115 real(kp), intent(in) :: trice
116 real(kp), intent(in) :: latvap
117 real(kp), intent(in) :: latice
118 real(kp), intent(in) :: rh2o
119 real(kp), intent(in) :: cpair
120 real(kp), intent(in) :: tmeltx
121!
122!---------------------------Local variables-----------------------------
123!
124 real(kp) t
125 integer n
126 integer lentbl
127 integer itype
128! 1 -> ice phase, no transition
129! -x -> ice phase, x degree transition
130 logical ip
131!
132!-----------------------------------------------------------------------
133!
134! Set es table parameters
135!
136 tmin = tmn
137 tmax = tmx
138 ttrice = trice
139 icephs = ip
140!
141! Set physical constants required for es calculation
142!
143 epsqs = epsil
144 hlatv = latvap
145 hlatf = latice
146 rgasv = rh2o
147 cp = cpair
148 tmelt = tmeltx
149!
150 lentbl = int(tmax-tmin+2.000001_kp)
151 if (lentbl .gt. plenest) then
152
153
154
155 write(*,*) "AHHH wv_sat"
156 stop
157 end if
158!
159! Begin building es table.
160! Check whether ice phase requested.
161! If so, set appropriate transition range for temperature
162!
163 if (icephs) then
164 if (ttrice /= 0.0_kp) then
165 itype = -ttrice
166 else
167 itype = 1
168 end if
169 else
170 itype = 0
171 end if
172!
173 t = tmin - 1.0_kp
174 do n=1,lentbl
175 t = t + 1.0_kp
176 call gffgch(t,estbl(n),tmelt,itype)
177 end do
178!
179 do n=lentbl+1,plenest
180 estbl(n) = -99999.0_kp
181 end do
182!
183! Table complete -- Set coefficients for polynomial approximation of
184! difference between saturation vapor press over water and saturation
185! pressure over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial
186! is valid in the range -40 < t < 0 (degrees C).
187!
188! --- Degree 5 approximation ---
189!
190 pcf(1) = 5.04469588506e-01_kp
191 pcf(2) = -5.47288442819e+00_kp
192 pcf(3) = -3.67471858735e-01_kp
193 pcf(4) = -8.95963532403e-03_kp
194 pcf(5) = -7.78053686625e-05_kp
195!
196! --- Degree 6 approximation ---
197!
198!-----pcf(1) = 7.63285250063e-02
199!-----pcf(2) = -5.86048427932e+00
200!-----pcf(3) = -4.38660831780e-01
201!-----pcf(4) = -1.37898276415e-02
202!-----pcf(5) = -2.14444472424e-04
203!-----pcf(6) = -1.36639103771e-06
204!
205
206!++jtb (comment out)
207! if (masterproc) then
208! !!write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'
209! end if
210!--jtb
211
212 return
213!
2149000 format('GESTBL: FATAL ERROR ********************************* &
215 &',/, ' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', &
216 & ' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, &
217 & ' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3)
218!
219 end subroutine gestbl
220
227 subroutine aqsat(t ,p ,es ,qs ,ii , ilen ,kk ,kstart ,kend )
228!------------------------------Arguments--------------------------------
229!
230! Input arguments
231!
232 integer, intent(in) :: ii
233 integer, intent(in) :: kk
234 integer, intent(in) :: ilen
235 integer, intent(in) :: kstart
236 integer, intent(in) :: kend
237 real(kp), intent(in) :: t(ii,kk)
238 real(kp), intent(in) :: p(ii,kk)
239!
240! Output arguments
241!
242 real(kp), intent(out) :: es(ii,kk)
243 real(kp), intent(out) :: qs(ii,kk)
244!
245!---------------------------Local workspace-----------------------------
246!
247 real(kp) omeps
248 integer i, k
249!
250!-----------------------------------------------------------------------
251!
252 omeps = 1.0_kp - epsqs
253 do k=kstart,kend
254 do i=1,ilen
255 es(i,k) = min(estblf(t(i,k)),p(i,k))
256!
257! Saturation specific humidity
258!
259 qs(i,k) = min(1.0_kp, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k)))
260!
261! The following check is to avoid the generation of negative values
262! that can occur in the upper stratosphere and mesosphere
263!
264! if (qs(i,k) < 0.0_kp) then
265! qs(i,k) = 1.0_kp
266! es(i,k) = p(i,k)
267! end if
268
269 end do
270 end do
271!
272 return
273 end subroutine aqsat
274
275!++xl
283 subroutine aqsat_water(t, p, es, qs, ii, ilen, kk, kstart,kend)
284!------------------------------Arguments--------------------------------
285!
286! Input arguments
287!
288 integer, intent(in) :: ii
289 integer, intent(in) :: kk
290 integer, intent(in) :: ilen
291 integer, intent(in) :: kstart
292 integer, intent(in) :: kend
293 real(kp), intent(in) :: t(ii,kk)
294 real(kp), intent(in) :: p(ii,kk)
295!
296! Output arguments
297!
298 real(kp), intent(out) :: es(ii,kk)
299 real(kp), intent(out) :: qs(ii,kk)
300!
301!---------------------------Local workspace-----------------------------
302!
303 real(kp) omeps
304 integer i, k
305!
306!-----------------------------------------------------------------------
307!
308 omeps = 1.0_kp - epsqs
309 do k=kstart,kend
310 do i=1,ilen
311! es(i,k) = estblf(t(i,k))
312#ifdef GEOS5
313 es(i,k) = min(polysvp(t(i,k),0), p(i,k))
314#endif
315#ifdef NEMS_GSM
316 es(i,k) = min(fpvsl(t(i,k)), p(i,k))
317#endif
318!
319! Saturation specific humidity
320!
321 qs(i,k) = min(1.0_kp, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k)))
322!
323! The following check is to avoid the generation of negative values
324! that can occur in the upper stratosphere and mesosphere
325!
326! if (qs(i,k) < 0.0_kp) then
327! qs(i,k) = 1.0_kp
328! es(i,k) = p(i,k)
329! end if
330 end do
331 end do
332!
333 return
334 end subroutine aqsat_water
335!--xl
336
348 subroutine aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend)
349!------------------------------Arguments--------------------------------
350!
351! Input arguments
352!
353 integer, intent(in) :: ii
354 integer, intent(in) :: ilen
355 integer, intent(in) :: kk
356 integer, intent(in) :: kstart
357 integer, intent(in) :: kend
358
359 real(kp), intent(in) :: t(ii,kk)
360 real(kp), intent(in) :: p(ii,kk)
361
362!
363! Output arguments
364!
365 real(kp), intent(out) :: es(ii,kk)
366 real(kp), intent(out) :: qs(ii,kk)
367 real(kp), intent(out) :: gam(ii,kk)
368!
369!---------------------------Local workspace-----------------------------
370!
371 logical lflg
372 integer i
373 integer k
374 real(kp) omeps
375 real(kp) trinv
376 real(kp) tc
377 real(kp) weight
378 real(kp) hltalt
379 real(kp) hlatsb
380 real(kp) hlatvp
381 real(kp) tterm
382 real(kp) desdt
383!
384!-----------------------------------------------------------------------
385!
386 omeps = 1.0_kp - epsqs
387 do k=kstart,kend
388 do i=1,ilen
389 es(i,k) = min(p(i,k), estblf(t(i,k)))
390!
391! Saturation specific humidity
392!
393 qs(i,k) = min(1.0_kp, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k)))
394!
395! The following check is to avoid the generation of negative qs
396! values which can occur in the upper stratosphere and mesosphere
397!
398!
399! if (qs(i,k) < 0.0_kp) then
400! qs(i,k) = 1.0_kp
401! es(i,k) = p(i,k)
402! end if
403 end do
404 end do
405!
406! "generalized" analytic expression for t derivative of es
407! accurate to within 1 percent for 173.16 < t < 373.16
408!
409 trinv = 0.0_kp
410 if ((.not. icephs) .or. (ttrice == 0.0_kp)) go to 10
411 trinv = 1.0_kp/ttrice
412!
413 do k=kstart,kend
414 do i=1,ilen
415!
416! Weighting of hlat accounts for transition from water to ice
417! polynomial expression approximates difference between es over
418! water and es over ice from 0 to -ttrice (C) (min of ttrice is
419! -40): required for accurate estimate of es derivative in transition
420! range from ice to water also accounting for change of hlatv with t
421! above freezing where constant slope is given by -2369 j/(kg c) =cpv - cw
422!
423 tc = t(i,k) - tmelt
424 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
425 weight = min(-tc*trinv,1.0_kp)
426 hlatsb = hlatv + weight*hlatf
427 hlatvp = hlatv - 2369.0_kp*tc
428 if (t(i,k) < tmelt) then
429 hltalt = hlatsb
430 else
431 hltalt = hlatvp
432 end if
433 if (lflg) then
434 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
435 & + tc*pcf(5))))
436 else
437 tterm = 0.0_kp
438 end if
439 desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) + tterm*trinv
440 gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) &
441 & - omeps*es(i,k)))
442 if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp
443 end do
444 end do
445!
446 go to 20
447!
448! No icephs or water to ice transition
449!
45010 do k=kstart,kend
451 do i=1,ilen
452!
453! Account for change of hlatv with t above freezing where
454! constant slope is given by -2369 j/(kg c) = cpv - cw
455!
456 hlatvp = hlatv - 2369.0_kp*(t(i,k)-tmelt)
457 if (icephs) then
458 hlatsb = hlatv + hlatf
459 else
460 hlatsb = hlatv
461 end if
462 if (t(i,k) < tmelt) then
463 hltalt = hlatsb
464 else
465 hltalt = hlatvp
466 end if
467 desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k))
468 gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) &
469 & - omeps*es(i,k)))
470 if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp
471 end do
472 end do
473!
47420 return
475 end subroutine aqsatd
476
484 subroutine vqsatd(t ,p ,es ,qs ,gam , len )
485!------------------------------Arguments--------------------------------
486!
487! Input arguments
488!
489 integer, intent(in) :: len
490 real(kp), intent(in) :: t(len)
491 real(kp), intent(in) :: p(len)
492!
493! Output arguments
494!
495 real(kp), intent(out) :: es(len)
496 real(kp), intent(out) :: qs(len)
497 real(kp), intent(out) :: gam(len)
498!
499!--------------------------Local Variables------------------------------
500!
501 logical lflg
502!
503 integer i
504!
505 real(kp) omeps
506 real(kp) trinv
507 real(kp) tc
508 real(kp) weight
509 real(kp) hltalt
510!
511 real(kp) hlatsb
512 real(kp) hlatvp
513 real(kp) tterm
514 real(kp) desdt
515!
516!-----------------------------------------------------------------------
517!
518 omeps = 1.0_kp - epsqs
519 do i=1,len
520 es(i) = min(estblf(t(i)), p(i))
521!
522! Saturation specific humidity
523!
524 qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
525!
526! The following check is to avoid the generation of negative
527! values that can occur in the upper stratosphere and mesosphere
528!
529 qs(i) = min(1.0_kp,qs(i))
530!
531! if (qs(i) < 0.0_kp) then
532! qs(i) = 1.0_kp
533! es(i) = p(i)
534! end if
535
536 end do
537!
538! "generalized" analytic expression for t derivative of es
539! accurate to within 1 percent for 173.16 < t < 373.16
540!
541 trinv = 0.0_kp
542 if ((.not. icephs) .or. (ttrice.eq.0.0_kp)) go to 10
543 trinv = 1.0_kp/ttrice
544 do i=1,len
545!
546! Weighting of hlat accounts for transition from water to ice
547! polynomial expression approximates difference between es over
548! water and es over ice from 0 to -ttrice (C) (min of ttrice is
549! -40): required for accurate estimate of es derivative in transition
550! range from ice to water also accounting for change of hlatv with t
551! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
552!
553 tc = t(i) - tmelt
554 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
555 weight = min(-tc*trinv,1.0_kp)
556 hlatsb = hlatv + weight*hlatf
557 hlatvp = hlatv - 2369.0_kp*tc
558 if (t(i) < tmelt) then
559 hltalt = hlatsb
560 else
561 hltalt = hlatvp
562 end if
563 if (lflg) then
564 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
565 & + tc*pcf(5))))
566 else
567 tterm = 0.0_kp
568 end if
569 desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
570 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
571 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
572 end do
573 return
574!
575! No icephs or water to ice transition
576!
57710 do i=1,len
578!
579! Account for change of hlatv with t above freezing where
580! constant slope is given by -2369 j/(kg c) = cpv - cw
581!
582 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
583 if (icephs) then
584 hlatsb = hlatv + hlatf
585 else
586 hlatsb = hlatv
587 end if
588 if (t(i) < tmelt) then
589 hltalt = hlatsb
590 else
591 hltalt = hlatvp
592 end if
593 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
594 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
595 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
596 end do
597!
598 return
599!
600 end subroutine vqsatd
601
602!++xl
605 subroutine vqsatd_water(t, p, es, qs, gam, len)
606
607!------------------------------Arguments--------------------------------
608!
609! Input arguments
610!
611 integer, intent(in) :: len
612 real(kp), intent(in) :: t(len)
613 real(kp), intent(in) :: p(len)
614
615!
616! Output arguments
617!
618 real(kp), intent(out) :: es(len)
619 real(kp), intent(out) :: qs(len)
620 real(kp), intent(out) :: gam(len)
621
622!
623!--------------------------Local Variables------------------------------
624!
625!
626 integer i
627!
628 real(kp) omeps
629 real(kp) hltalt
630!
631 real(kp) hlatsb
632 real(kp) hlatvp
633 real(kp) desdt
634!
635!-----------------------------------------------------------------------
636!
637 omeps = 1.0_kp - epsqs
638 do i=1,len
639#ifdef NEMS_GSM
640 es(i) = min(fpvsl(t(i)), p(i))
641#else
642 es(i) = min(polysvp(t(i),0), p(i))
643#endif
644!
645! Saturation specific humidity
646!
647 qs(i) = min(1.0_kp, epsqs*es(i) / (p(i)-omeps*es(i)))
648!
649! The following check is to avoid the generation of negative
650! values that can occur in the upper stratosphere and mesosphere
651!
652! qs(i) = min(1.0_kp,qs(i))
653!
654! if (qs(i) < 0.0_kp) then
655! qs(i) = 1.0_kp
656! es(i) = p(i)
657! end if
658
659 end do
660!
661! No icephs or water to ice transition
662!
663 do i=1,len
664!
665! Account for change of hlatv with t above freezing where
666! constant slope is given by -2369 j/(kg c) = cpv - cw
667!
668 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
669 hlatsb = hlatv
670 if (t(i) < tmelt) then
671 hltalt = hlatsb
672 else
673 hltalt = hlatvp
674 end if
675 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
676 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
677 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
678 end do
679!
680 return
681!
682 end subroutine vqsatd_water
683
685 function polysvp (T,typ)
686
691
692!!DONIFF Changed to Murphy and Koop (2005) (03/04/14)
693
694
695 real(kp) dum
696
697 real(kp) t,polysvp
698
699 integer typ
700
701
702 if (.true.) then
703!ice
704 if (typ == 1) then
705 polysvp = murphykoop_svp_ice(t)
706 end if
707 if (typ == 0) then
709 end if
710
711 else
712
713! ice
714 if (typ.eq.1) then
715
716
717
718 polysvp = 10._kp**(-9.09718_kp*(273.16_kp/t-1._kp)-3.56654_kp* &
719 & log10(273.16_kp/t)+0.876793_kp*(1._kp-t/273.16_kp)+ &
720 & log10(6.1071_kp))*100._kp
721
722 end if
723
724
725
726 if (typ.eq.0) then
727 polysvp = 10._kp**(-7.90298_kp*(373.16_kp/t-1._kp)+ 5.02808_kp* &
728 &log10(373.16_kp/t)- 1.3816e-7_kp*(10._kp**(11.344_kp*(1._kp-t/ &
729 &373.16_kp))-1._kp)+ 8.1328e-3_kp*(10._kp**(-3.49149_kp*(373.16_kp/ &
730 &t-1._kp))-1._kp)+ log10(1013.246_kp))*100._kp
731 end if
732
733 end if
734
735 end function polysvp
736!--xl
737
738
739
740 integer function fqsatd(t ,p ,es ,qs ,gam , len )
741
742
743
744
745
746 integer, intent(in) :: len
747 real(kp), intent(in) :: t(len)
748 real(kp), intent(in) :: p(len)
749
750 real(kp), intent(out) :: es(len)
751 real(kp), intent(out) :: qs(len)
752 real(kp), intent(out) :: gam(len)
753
754 call vqsatd(t ,p ,es ,qs ,gam , len )
755 fqsatd = 1
756 return
757 end function fqsatd
758
759 real(kp) function qsat_water(t,p)
760
761 real(kp) t
762 real(kp) p
763 real(kp) es
764 real(kp) ps, ts, e1, e2, f1, f2, f3, f4, f5, f
765
766
767
768
769
770 ps = 1013.246_kp
771 ts = 373.16_kp
772 e1 = 11.344_kp*(1.0_kp - t/ts)
773 e2 = -3.49149_kp*(ts/t - 1.0_kp)
774 f1 = -7.90298_kp*(ts/t - 1.0_kp)
775 f2 = 5.02808_kp*log10(ts/t)
776 f3 = -1.3816_kp*(10.0_kp**e1 - 1.0_kp)/10000000.0_kp
777 f4 = 8.1328_kp*(10.0_kp**e2 - 1.0_kp)/1000.0_kp
778 f5 = log10(ps)
779 f = f1 + f2 + f3 + f4 + f5
780 es = (10.0_kp**f)*100.0_kp
781
782 qsat_water = epsqs*es/(p-(1.-epsqs)*es)
783 if(qsat_water < 0.) qsat_water = 1.
784
785 end function qsat_water
786
789 subroutine vqsat_water(t,p,qsat_water,len)
790
791 integer, intent(in) :: len
792 real(kp) t(len)
793 real(kp) p(len)
794 real(kp) qsat_water(len)
795 real(kp) es
796 real(kp), parameter :: t0inv = 1._kp/273._kp
797 real(kp) coef
798 integer :: i
799
800 coef = hlatv/rgasv
801 do i=1,len
802 es = 611._kp*exp(coef*(t0inv-1./t(i)))
803 qsat_water(i) = epsqs*es/(p(i)-(1.-epsqs)*es)
804 if(qsat_water(i) < 0.) qsat_water(i) = 1.
805 enddo
806
807 return
808
809 end subroutine vqsat_water
810
812 real(kp) function qsat_ice(t,p)
813
814 real(kp) t
815 real(kp) p
816 real(kp) es
817 real(kp), parameter :: t0inv = 1._kp/273._kp
818 es = 611.*exp((hlatv+hlatf)/rgasv*(t0inv-1./t))
819 qsat_ice = epsqs*es/(p-(1.-epsqs)*es)
820 if(qsat_ice < 0.) qsat_ice = 1.
821
822 end function qsat_ice
823
825 subroutine vqsat_ice(t,p,qsat_ice,len)
826
827 integer,intent(in) :: len
828 real(kp) t(len)
829 real(kp) p(len)
830 real(kp) qsat_ice(len)
831 real(kp) es
832 real(kp), parameter :: t0inv = 1._kp/273._kp
833 real(kp) coef
834 integer :: i
835
836 coef = (hlatv+hlatf)/rgasv
837 do i=1,len
838 es = 611.*exp(coef*(t0inv-1./t(i)))
839 qsat_ice(i) = epsqs*es/(p(i)-(1.-epsqs)*es)
840 if(qsat_ice(i) < 0.) qsat_ice(i) = 1.
841 enddo
842
843 return
844
845 end subroutine vqsat_ice
846
847! Sungsu
848! Below two subroutines (vqsatd2_water,vqsatd2_water_single) are by Sungsu
849! Replace 'gam -> dqsdt'
850! Sungsu
851
853 subroutine vqsatd2_water(t ,p ,es ,qs ,dqsdt , len )
854
855!------------------------------Arguments--------------------------------
856!
857! Input arguments
858!
859 integer, intent(in) :: len
860 real(kp), intent(in) :: t(len)
861 real(kp), intent(in) :: p(len)
862
863!
864! Output arguments
865!
866 real(kp), intent(out) :: es(len)
867 real(kp), intent(out) :: qs(len)
868
869
870 real(kp), intent(out) :: dqsdt(len)
871
872
873!
874!--------------------------Local Variables------------------------------
875!
876!
877 integer i
878!
879 real(kp) omeps
880 real(kp) hltalt
881!
882 real(kp) hlatsb
883 real(kp) hlatvp
884 real(kp) desdt
885
886
887 real(kp) gam(len)
888
889
890!
891!-----------------------------------------------------------------------
892!
893 omeps = 1.0_kp - epsqs
894 do i=1,len
895#ifdef GEOS5
896 es(i) = min(polysvp(t(i),0), p(i))
897#endif
898#ifdef NEMS_GSM
899 es(i) = min(fpvsl(t(i)), p(i))
900#endif
901!
902! Saturation specific humidity
903!
904 qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
905!
906! The following check is to avoid the generation of negative
907! values that can occur in the upper stratosphere and mesosphere
908!
909 qs(i) = min(1.0_kp,qs(i))
910!
911! if (qs(i) < 0.0_kp) then
912! qs(i) = 1.0_kp
913! es(i) = p(i)
914! end if
915
916 end do
917!
918! No icephs or water to ice transition
919!
920 do i=1,len
921!
922! Account for change of hlatv with t above freezing where
923! constant slope is given by -2369 j/(kg c) = cpv - cw
924!
925 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
926 hlatsb = hlatv
927 if (t(i) < tmelt) then
928 hltalt = hlatsb
929 else
930 hltalt = hlatvp
931 end if
932 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
933 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
934 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
935
936 dqsdt(i) = (cp/hltalt)*gam(i)
937
938 end do
939!
940 return
941!
942 end subroutine vqsatd2_water
943
945 subroutine vqsatd2_water_single(t ,p ,es ,qs ,dqsdt)
946
947!------------------------------Arguments--------------------------------
948!
949! Input arguments
950!
951
952 real(kp), intent(in) :: t, p
953
954!
955! Output arguments
956!
957 real(kp), intent(out) :: es, qs, dqsdt
958!
959!--------------------------Local Variables------------------------------
960!
961! integer i
962!
963 real(kp) omeps, hltalt, hlatsb, hlatvp, desdt, gam
964!
965!-----------------------------------------------------------------------
966!
967 omeps = 1.0_kp - epsqs
968! do i=1,len
969#ifdef GEOS5
970 es = min(p, polysvp(t,0))
971#endif
972#ifdef NEMS_GSM
973 es = min(p, fpvsl(t))
974#endif
975!
976! Saturation specific humidity
977!
978 qs = min(1.0_kp, epsqs*es/(p-omeps*es))
979!
980! The following check is to avoid the generation of negative
981! values that can occur in the upper stratosphere and mesosphere
982!
983! if (qs < 0.0_kp) then
984! qs = 1.0_kp
985! es = p
986! end if
987! end do
988!
989! No icephs or water to ice transition
990!
991! do i=1,len
992!
993! Account for change of hlatv with t above freezing where
994! constant slope is given by -2369 j/(kg c) = cpv - cw
995!
996 hlatvp = hlatv - 2369.0_kp*(t-tmelt)
997 hlatsb = hlatv
998 if (t < tmelt) then
999 hltalt = hlatsb
1000 else
1001 hltalt = hlatvp
1002 end if
1003 desdt = hltalt*es/(rgasv*t*t)
1004 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1005 if (qs >= 1.0_kp) gam = 0.0_kp
1006
1007 dqsdt = (cp/hltalt)*gam
1008
1009! end do
1010!
1011 return
1012!
1013 end subroutine vqsatd2_water_single
1014
1016 subroutine vqsatd2(t ,p ,es ,qs ,dqsdt , len )
1017!-----------------------------------------------------------------------
1018! Sungsu : This is directly copied from 'vqsatd' but 'dqsdt' is output
1019! instead of gam for use in Sungsu's equilibrium stratiform
1020! macrophysics scheme.
1021!
1022! Purpose:
1023! Utility procedure to look up and return saturation vapor pressure from
1024! precomputed table, calculate and return saturation specific humidity
1025! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT). The same
1026! function as qsatd, but operates on vectors of temperature and pressure
1027!
1028! Method:
1029!
1030! Author: J. Hack
1031!
1032!------------------------------Arguments--------------------------------
1033!
1034! Input arguments
1035!
1036 integer, intent(in) :: len
1037 real(kp), intent(in) :: t(len)
1038 real(kp), intent(in) :: p(len)
1039!
1040! Output arguments
1041!
1042 real(kp), intent(out) :: es(len)
1043 real(kp), intent(out) :: qs(len)
1044
1045
1046 real(kp), intent(out) :: dqsdt(len)
1047
1048
1049!
1050!--------------------------Local Variables------------------------------
1051!
1052 logical lflg
1053!
1054 integer i
1055!
1056 real(kp) omeps
1057 real(kp) trinv
1058 real(kp) tc
1059 real(kp) weight
1060 real(kp) hltalt
1061!
1062 real(kp) hlatsb
1063 real(kp) hlatvp
1064 real(kp) tterm
1065 real(kp) desdt
1066
1067
1068 real(kp) gam(len)
1069
1070!
1071!-----------------------------------------------------------------------
1072!
1073 omeps = 1.0_kp - epsqs
1074 do i=1,len
1075#ifdef GEOS5
1076 es(i) = min(p(i), estblf(t(i)))
1077#endif
1078#ifdef NEMS_GSM
1079 es(i) = min(p(i), fpvsi(t(i)))
1080#endif
1081!
1082! Saturation specific humidity
1083!
1084 qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
1085!
1086! The following check is to avoid the generation of negative
1087! values that can occur in the upper stratosphere and mesosphere
1088!
1089 qs(i) = min(1.0_kp,qs(i))
1090!
1091! if (qs(i) < 0.0_kp) then
1092! qs(i) = 1.0_kp
1093! es(i) = p(i)
1094! end if
1095 end do
1096!
1097! "generalized" analytic expression for t derivative of es
1098! accurate to within 1 percent for 173.16 < t < 373.16
1099!
1100 trinv = 0.0_kp
1101 if ((.not. icephs) .or. (ttrice == 0.0_kp)) go to 10
1102 trinv = 1.0_kp/ttrice
1103 do i=1,len
1104!
1105! Weighting of hlat accounts for transition from water to ice
1106! polynomial expression approximates difference between es over
1107! water and es over ice from 0 to -ttrice (C) (min of ttrice is
1108! -40): required for accurate estimate of es derivative in transition
1109! range from ice to water also accounting for change of hlatv with t
1110! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
1111!
1112 tc = t(i) - tmelt
1113 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
1114 weight = min(-tc*trinv,1.0_kp)
1115 hlatsb = hlatv + weight*hlatf
1116 hlatvp = hlatv - 2369.0_kp*tc
1117 if (t(i) < tmelt) then
1118 hltalt = hlatsb
1119 else
1120 hltalt = hlatvp
1121 end if
1122 if (lflg) then
1123 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
1124 & + tc*pcf(5))))
1125 else
1126 tterm = 0.0_kp
1127 end if
1128 desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
1129 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
1130 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
1131
1132 dqsdt(i) = (cp/hltalt)*gam(i)
1133
1134 end do
1135 return
1136!
1137! No icephs or water to ice transition
1138!
113910 do i=1,len
1140!
1141! Account for change of hlatv with t above freezing where
1142! constant slope is given by -2369 j/(kg c) = cpv - cw
1143!
1144 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
1145 if (icephs) then
1146 hlatsb = hlatv + hlatf
1147 else
1148 hlatsb = hlatv
1149 end if
1150 if (t(i) < tmelt) then
1151 hltalt = hlatsb
1152 else
1153 hltalt = hlatvp
1154 end if
1155 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
1156 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
1157 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
1158
1159 dqsdt(i) = (cp/hltalt)*gam(i)
1160
1161 end do
1162!
1163 return
1164!
1165 end subroutine vqsatd2
1166
1167
1168! Below routine is by Sungsu
1169
1171 subroutine vqsatd2_single(t ,p ,es ,qs ,dqsdt)
1172!-----------------------------------------------------------------------
1173! Sungsu : This is directly copied from 'vqsatd' but 'dqsdt' is output
1174! instead of gam for use in Sungsu's equilibrium stratiform
1175! macrophysics scheme.
1176!
1177! Purpose:
1178! Utility procedure to look up and return saturation vapor pressure from
1179! precomputed table, calculate and return saturation specific humidity
1180! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT). The same
1181! function as qsatd, but operates on vectors of temperature and pressure
1182!
1183! Method:
1184!
1185! Author: J. Hack
1186!
1187!------------------------------Arguments--------------------------------
1188!
1189! Input arguments
1190!
1191 real(kp), intent(in) :: t, p
1192!
1193! Output arguments
1194!
1195 real(kp), intent(out) :: es, qs, dqsdt
1196!
1197!--------------------------Local Variables------------------------------
1198!
1199 logical lflg
1200!
1201! integer i ! index for vector calculations
1202!
1203 real(kp) omeps
1204 real(kp) trinv
1205 real(kp) tc
1206 real(kp) weight
1207 real(kp) hltalt
1208!
1209 real(kp) hlatsb
1210 real(kp) hlatvp
1211 real(kp) tterm
1212 real(kp) desdt
1213
1214 real(kp) gam
1215
1216!
1217!-----------------------------------------------------------------------
1218!
1219 omeps = 1.0_kp - epsqs
1220
1221! do i=1,len
1222
1223#ifdef GEOS5
1224 es = estblf(t)
1225#endif
1226#ifdef NEMS_GSM
1227 es = min(fpvs(t), p)
1228#endif
1229!
1230! Saturation specific humidity
1231!
1232 qs = epsqs*es/(p - omeps*es)
1233!
1234! The following check is to avoid the generation of negative
1235! values that can occur in the upper stratosphere and mesosphere
1236!
1237 qs = min(1.0_kp,qs)
1238!
1239! if (qs < 0.0_kp) then
1240! qs = 1.0_kp
1241! es = p
1242! end if
1243
1244! end do
1245!
1246! "generalized" analytic expression for t derivative of es
1247! accurate to within 1 percent for 173.16 < t < 373.16
1248!
1249 trinv = 0.0_kp
1250 if ((.not. icephs) .or. (ttrice == 0.0_kp)) go to 10
1251 trinv = 1.0_kp/ttrice
1252
1253! do i=1,len
1254!
1255! Weighting of hlat accounts for transition from water to ice
1256! polynomial expression approximates difference between es over
1257! water and es over ice from 0 to -ttrice (C) (min of ttrice is
1258! -40): required for accurate estimate of es derivative in transition
1259! range from ice to water also accounting for change of hlatv with t
1260! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
1261!
1262 tc = t - tmelt
1263 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
1264 weight = min(-tc*trinv,1.0_kp)
1265 hlatsb = hlatv + weight*hlatf
1266 hlatvp = hlatv - 2369.0_kp*tc
1267 if (t < tmelt) then
1268 hltalt = hlatsb
1269 else
1270 hltalt = hlatvp
1271 end if
1272 if (lflg) then
1273 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
1274 & + tc*pcf(5))))
1275 else
1276 tterm = 0.0_kp
1277 end if
1278 desdt = hltalt*es/(rgasv*t*t) + tterm*trinv
1279 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1280 if (qs == 1.0_kp) gam = 0.0_kp
1281
1282 dqsdt = (cp/hltalt)*gam
1283
1284! end do
1285 return
1286!
1287! No icephs or water to ice transition
1288!
1289
129010 continue
1291
1292!10 do i=1,len
1293!
1294! Account for change of hlatv with t above freezing where
1295! constant slope is given by -2369 j/(kg c) = cpv - cw
1296!
1297 hlatvp = hlatv - 2369.0_kp*(t-tmelt)
1298 if (icephs) then
1299 hlatsb = hlatv + hlatf
1300 else
1301 hlatsb = hlatv
1302 end if
1303 if (t < tmelt) then
1304 hltalt = hlatsb
1305 else
1306 hltalt = hlatvp
1307 end if
1308 desdt = hltalt*es/(rgasv*t*t)
1309 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1310 if (qs == 1.0_kp) gam = 0.0_kp
1311
1312 dqsdt = (cp/hltalt)*gam
1313
1314
1315! end do
1316!
1317 return
1318!
1319 end subroutine vqsatd2_single
1320
1321!----------------------------------------------------------------------
1322
1323!----------------------------------------------------------------------
1324
1326 subroutine gffgch(t ,es ,tmelt ,itype )
1327!-----------------------------------------------------------------------
1328!
1329! Purpose:
1330! Computes saturation vapor pressure over water and/or over ice using
1331! Goff & Gratch (1946) relationships.
1332! <Say what the routine does>
1333!
1334! Method:
1335! T (temperature), and itype are input parameters, while es (saturation
1336! vapor pressure) is an output parameter. The input parameter itype
1337! serves two purposes: a value of zero indicates that saturation vapor
1338! pressures over water are to be returned (regardless of temperature),
1339! while a value of one indicates that saturation vapor pressures over
1340! ice should be returned when t is less than freezing degrees. If itype
1341! is negative, its absolute value is interpreted to define a temperature
1342! transition region below freezing in which the returned
1343! saturation vapor pressure is a weighted average of the respective ice
1344! and water value. That is, in the temperature range 0 => -itype
1345! degrees c, the saturation vapor pressures are assumed to be a weighted
1346! average of the vapor pressure over supercooled water and ice (all
1347! water at 0 c; all ice at -itype c). Maximum transition range => 40 c
1348!
1349! Author: J. Hack
1350!
1351!-----------------------------------------------------------------------
1352
1353
1354
1355
1356
1357 implicit none
1358!------------------------------Arguments--------------------------------
1359!
1360! Input arguments
1361!
1362 real(kp), intent(in) :: t ,tmelt
1363!
1364! Output arguments
1365!
1366 integer, intent(inout) :: itype
1367
1368 real(kp), intent(out) :: es
1369!
1370!---------------------------Local variables-----------------------------
1371!
1372 real(kp) e1
1373 real(kp) e2
1374 real(kp) eswtr
1375 real(kp) f
1376 real(kp) f1
1377 real(kp) f2
1378 real(kp) f3
1379 real(kp) f4
1380 real(kp) f5
1381 real(kp) ps
1382 real(kp) t0
1383 real(kp) term1
1384 real(kp) term2
1385 real(kp) term3
1386 real(kp) tr
1387 real(kp) ts
1388 real(kp) weight
1389 integer itypo
1390!
1391!-----------------------------------------------------------------------
1392!
1393! Check on whether there is to be a transition region for es
1394!
1395 if (itype < 0) then
1396 tr = abs(real(itype,kp))
1397 itypo = itype
1398 itype = 1
1399 else
1400 tr = 0.0_kp
1401 itypo = itype
1402 end if
1403 if (tr > 40.0_kp) then
1404 write(iulog,900) tr
1405
1406 end if
1407!
1408 if(t < (tmelt - tr) .and. itype == 1) go to 10
1409!
1410! Water
1411!
1412 ps = 1013.246_kp
1413 ts = 373.16_kp
1414 e1 = 11.344_kp*(1.0_kp - t/ts)
1415 e2 = -3.49149_kp*(ts/t - 1.0_kp)
1416 f1 = -7.90298_kp*(ts/t - 1.0_kp)
1417 f2 = 5.02808_kp*log10(ts/t)
1418 f3 = -1.3816_kp*(10.0_kp**e1 - 1.0_kp)/10000000.0_kp
1419 f4 = 8.1328_kp*(10.0_kp**e2 - 1.0_kp)/1000.0_kp
1420 f5 = log10(ps)
1421 f = f1 + f2 + f3 + f4 + f5
1422 es = (10.0_kp**f)*100.0_kp
1423 eswtr = es
1424!
1425 if(t >= tmelt .or. itype == 0) go to 20
1426!
1427! Ice
1428!
142910 continue
1430 t0 = tmelt
1431 term1 = 2.01889049_kp/(t0/t)
1432 term2 = 3.56654_kp*log(t0/t)
1433 term3 = 20.947031_kp*(t0/t)
1434 es = 575.185606e10_kp*exp(-(term1 + term2 + term3))
1435!
1436 if (t < (tmelt - tr)) go to 20
1437!
1438! Weighted transition between water and ice
1439!
1440 weight = min((tmelt - t)/tr,1.0_kp)
1441 es = weight*es + (1.0_kp - weight)*eswtr
1442!
144320 continue
1444 itype = itypo
1445 return
1446!
1447900 format('GFFGCH: FATAL ERROR ******************************',/, &
1448 & 'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR', ' PRESSURE, &
1449 & TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF', ' 40.0 DEGREES C',/, &
1450 & ' TR = ',f7.2)
1451!
1452 end subroutine gffgch
1453
1456 function murphykoop_svp_water(tx) result(es)
1457 real(kp), intent(in) :: tx
1458 real(kp) :: es
1459 real(kp):: t
1460
1461 t=min(tx, 332.0_kp)
1462 t=max(123.0_kp, tx)
1463
1464 es = exp(54.842763_kp - (6763.22_kp / t) - (4.210_kp * log(t)) + &
1465 & (0.000367_kp * t) + (tanh(0.0415_kp * (t - 218.8_kp)) * &
1466 & (53.878_kp - (1331.22_kp / t) - (9.44523_kp * log(t)) + &
1467 & 0.014025_kp * t)))
1468
1469 end function murphykoop_svp_water
1470
1471 function murphykoop_svp_ice(tx) result(es)
1472 real(kp), intent(in) :: tx
1473 real(kp) :: t
1474 real(kp) :: es
1475
1476 t=max(100.0_kp, tx)
1477 t=min(274.0_kp, tx)
1478
1479
1480 es = exp(9.550426_kp - (5723.265_kp / t) + (3.53068_kp * &
1481 & log(t)) - (0.00728332_kp * t))
1482
1483 end function murphykoop_svp_ice
1484!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1486 subroutine vqsatd2_ice_single(t ,p ,es ,qs ,dqsdt)
1487
1488!------------------------------Arguments--------------------------------
1489!
1490! Input arguments
1491!
1492 real(kp), intent(in) :: t, p
1493!
1494! Output arguments
1495!
1496 real(kp), intent(out) :: es, qs, dqsdt
1497
1498!
1499!--------------------------Local Variables------------------------------
1500!
1501! integer i
1502!
1503 real(kp) omeps, hltalt, hlatsb, hlatvp, desdt, gam
1504!
1505!-----------------------------------------------------------------------
1506!
1507 omeps = 1.0_kp - epsqs
1508! do i=1,len
1509#ifdef GEOS5
1510 es = min(polysvp(t,1),p)
1511#endif
1512#ifdef NEMS_GSM
1513 es = min(fpvsi(t),p)
1514#endif
1515!
1516! Saturation specific humidity
1517!
1518 qs = min(1.0_kp, epsqs*es/(p-omeps*es))
1519!
1520! The following check is to avoid the generation of negative
1521! values that can occur in the upper stratosphere and mesosphere
1522!
1523! if (qs < 0.0_kp) then
1524! qs = 1.0_kp
1525! es = p
1526! end if
1527! end do
1528!
1529! No icephs or water to ice transition
1530!
1531! do i=1,len
1532!
1533! Account for change of hlatv with t above freezing where
1534! constant slope is given by -2369 j/(kg c) = cpv - cw
1535!
1536 hltalt = hlatv + hlatf
1537 desdt = hltalt*es/(rgasv*t*t)
1538 if (qs < 1.0_kp) then
1539 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1540 else
1541 gam = 0.0_kp
1542 endif
1543
1544 dqsdt = (cp/hltalt)*gam
1545
1546! end do
1547!
1548 return
1549!
1550 end subroutine vqsatd2_ice_single
1551
1552!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1553
1554 end module wv_saturation
1555!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine, public aqsat(t, p, es, qs, ii, ilen, kk, kstart, kend)
This subroutine is the utility procedure to look up and returen saturation vapor pressure from precom...
subroutine, public aqsat_water(t, p, es, qs, ii, ilen, kk, kstart, kend)
This subroutine includes the utility procedure to look up and return saturation vapor pressure from p...
subroutine, public aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend)
This subroutine include the utility procedure to look up and returen saturation vapor pressure from p...
subroutine, public vqsatd(t, p, es, qs, gam, len)
This subroutine is the utility procedure to look up and return saturation vapor pressure from precomp...
subroutine, public vqsatd_water(t, p, es, qs, gam, len)
subroutine, public vqsat_water(t, p, qsat_water, len)
This subroutine.
subroutine, public vqsatd2_water_single(t, p, es, qs, dqsdt)
real(kp) function, public qsat_ice(t, p)
real(kp) function murphykoop_svp_water(tx)
DONIF USe Murphy and Koop (2005) (Written by Andrew Gettelman)
subroutine, public vqsatd2(t, p, es, qs, dqsdt, len)
real(kp) function, public polysvp(t, typ)
subroutine gffgch(t, es, tmelt, itype)
subroutine, public vqsatd2_ice_single(t, p, es, qs, dqsdt)
subroutine, public gestbl(tmn, tmx, trice, ip, epsil, latvap, latice, rh2o, cpair, tmeltx)
This subroutine builds saturation vapor pressure table for later procedure.
subroutine, public vqsatd2_water(t, p, es, qs, dqsdt, len)
subroutine, public vqsatd2_single(t, p, es, qs, dqsdt)
subroutine, public vqsat_ice(t, p, qsat_ice, len)
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26