UFO
ufo_utils_mod.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! (C) British Crown Copyright 2020 Met Office
3 !
4 ! This software is licensed under the terms of the Apache Licence Version 2.0
5 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 !-------------------------------------------------------------------------------
7 
8 !> Fortran module with various useful routines
9 
11 
12 use fckit_log_module, only : fckit_log
13 use kinds, only: kind_real
14 use missing_values_mod, only: missing_value
15 use ufo_constants_mod, only : zero, half, one, two, &
17 
18 implicit none
19 private
20 
21 public ops_qsat
22 public ops_qsatwat
23 public ops_satrad_qsplit
24 public ops_cholesky
26 public invertmatrix
27 public cmp_strings
28 
29 contains
30 
31 !-------------------------------------------------------------------------------
32 !> Calculate the Saturation Specific Humidity Scheme (Qsat): Vapour to Liquid/Ice.
33 !!
34 !! \details Heritage: Ops_Qsat.inc
35 !!
36 !! Returns a saturation mixing ratio given a temperature and pressure
37 !! using saturation vapour pressures caluclated using the Goff-Gratch
38 !! formulae, adopted by the WMO as taken from Landolt-Bornstein, 1987
39 !! Numerical data and Functional relationships in Science and
40 !! Technology. Group V/Vol 4B Meteorology. Physical and Chemical
41 !! properties of Air, P35.
42 !!
43 !! Value in the lookup table are over water above 0 degrees C and over
44 !! ice below this temperatures.
45 !!
46 !! Method: <br>
47 !! uses lookup tables to find eSAT, calculates qSAT directly from that.
48 !!
49 !! \author Met Office
50 !!
51 !! \date 09/06/2020: Created
52 !!
53 subroutine ops_qsat (QS, &
54  T, &
55  P, &
56  npnts)
57 
58 implicit none
59 
60 ! subroutine arguments:
61 integer, intent(in) :: npnts !< Points being processed by qSAT scheme.
62 real(kind=kind_real), intent(in) :: t(npnts) !< Temperature (K)
63 real(kind=kind_real), intent(in) :: p(npnts) !< Pressure (Pa).
64 real(kind=kind_real), intent(inout) :: qs(npnts) !< Saturation mixing ratio (KG/KG)
65 
66 ! Local declarations:
67 real(kind=kind_real), parameter :: one_minus_epsilon = one - epsilon
68 real(kind=kind_real), parameter :: t_low = 183.15_kind_real ! Lowest temperature for which look-up table is valid
69 real(kind=kind_real), parameter :: t_high = 338.15_kind_real ! Highest temperature for which look-up table is valid
70 real(kind=kind_real), parameter :: delta_t = 0.1_kind_real ! Temperature increment of look-up table
71 integer, parameter :: n = ((t_high - t_low + (delta_t * half)) / delta_t) + one ! Size of lookup-table (gives 1551)
72 integer :: itable
73 real(kind=kind_real) :: atable
74 real(kind=kind_real) :: fsubw ! Converts from sat vapour pressure in pure water to pressure in air
75 real(kind=kind_real) :: tt
76 integer :: i
77 integer :: ies
78 real(kind=kind_real) :: es(0:n + 1) ! Table of saturation water vapour pressure (PA)
79 character(len=*), parameter :: routinename = "Ops_Qsat"
80 
81 ! Note: 0 element is a repeat of 1st element to cater for special case
82 ! of low temperatures (.LE.T_low) for which the array index is
83 ! rounded down due to machine precision.
84 
85 data (es(ies), ies= 0, 95) / 0.966483d-02, &
86 0.966483d-02,0.984279d-02,0.100240d-01,0.102082d-01,0.103957d-01, &
87 0.105865d-01,0.107803d-01,0.109777d-01,0.111784d-01,0.113825d-01, &
88 0.115902d-01,0.118016d-01,0.120164d-01,0.122348d-01,0.124572d-01, &
89 0.126831d-01,0.129132d-01,0.131470d-01,0.133846d-01,0.136264d-01, &
90 0.138724d-01,0.141225d-01,0.143771d-01,0.146356d-01,0.148985d-01, &
91 0.151661d-01,0.154379d-01,0.157145d-01,0.159958d-01,0.162817d-01, &
92 0.165725d-01,0.168680d-01,0.171684d-01,0.174742d-01,0.177847d-01, &
93 0.181008d-01,0.184216d-01,0.187481d-01,0.190801d-01,0.194175d-01, &
94 0.197608d-01,0.201094d-01,0.204637d-01,0.208242d-01,0.211906d-01, &
95 0.215631d-01,0.219416d-01,0.223263d-01,0.227172d-01,0.231146d-01, &
96 0.235188d-01,0.239296d-01,0.243465d-01,0.247708d-01,0.252019d-01, &
97 0.256405d-01,0.260857d-01,0.265385d-01,0.269979d-01,0.274656d-01, &
98 0.279405d-01,0.284232d-01,0.289142d-01,0.294124d-01,0.299192d-01, &
99 0.304341d-01,0.309571d-01,0.314886d-01,0.320285d-01,0.325769d-01, &
100 0.331348d-01,0.337014d-01,0.342771d-01,0.348618d-01,0.354557d-01, &
101 0.360598d-01,0.366727d-01,0.372958d-01,0.379289d-01,0.385717d-01, &
102 0.392248d-01,0.398889d-01,0.405633d-01,0.412474d-01,0.419430d-01, &
103 0.426505d-01,0.433678d-01,0.440974d-01,0.448374d-01,0.455896d-01, &
104 0.463545d-01,0.471303d-01,0.479191d-01,0.487190d-01,0.495322d-01/
105  data (es(ies),ies= 96,190) / &
106 0.503591d-01,0.511977d-01,0.520490d-01,0.529145d-01,0.537931d-01, &
107 0.546854d-01,0.555924d-01,0.565119d-01,0.574467d-01,0.583959d-01, &
108 0.593592d-01,0.603387d-01,0.613316d-01,0.623409d-01,0.633655d-01, &
109 0.644053d-01,0.654624d-01,0.665358d-01,0.676233d-01,0.687302d-01, &
110 0.698524d-01,0.709929d-01,0.721490d-01,0.733238d-01,0.745180d-01, &
111 0.757281d-01,0.769578d-01,0.782061d-01,0.794728d-01,0.807583d-01, &
112 0.820647d-01,0.833905d-01,0.847358d-01,0.861028d-01,0.874882d-01, &
113 0.888957d-01,0.903243d-01,0.917736d-01,0.932464d-01,0.947407d-01, &
114 0.962571d-01,0.977955d-01,0.993584d-01,0.100942d+00,0.102551d+00, &
115 0.104186d+00,0.105842d+00,0.107524d+00,0.109231d+00,0.110963d+00, &
116 0.112722d+00,0.114506d+00,0.116317d+00,0.118153d+00,0.120019d+00, &
117 0.121911d+00,0.123831d+00,0.125778d+00,0.127755d+00,0.129761d+00, &
118 0.131796d+00,0.133863d+00,0.135956d+00,0.138082d+00,0.140241d+00, &
119 0.142428d+00,0.144649d+00,0.146902d+00,0.149190d+00,0.151506d+00, &
120 0.153859d+00,0.156245d+00,0.158669d+00,0.161126d+00,0.163618d+00, &
121 0.166145d+00,0.168711d+00,0.171313d+00,0.173951d+00,0.176626d+00, &
122 0.179342d+00,0.182096d+00,0.184893d+00,0.187724d+00,0.190600d+00, &
123 0.193518d+00,0.196473d+00,0.199474d+00,0.202516d+00,0.205604d+00, &
124 0.208730d+00,0.211905d+00,0.215127d+00,0.218389d+00,0.221701d+00/
125  data (es(ies),ies=191,285) / &
126 0.225063d+00,0.228466d+00,0.231920d+00,0.235421d+00,0.238976d+00, &
127 0.242580d+00,0.246232d+00,0.249933d+00,0.253691d+00,0.257499d+00, &
128 0.261359d+00,0.265278d+00,0.269249d+00,0.273274d+00,0.277358d+00, &
129 0.281498d+00,0.285694d+00,0.289952d+00,0.294268d+00,0.298641d+00, &
130 0.303078d+00,0.307577d+00,0.312135d+00,0.316753d+00,0.321440d+00, &
131 0.326196d+00,0.331009d+00,0.335893d+00,0.340842d+00,0.345863d+00, &
132 0.350951d+00,0.356106d+00,0.361337d+00,0.366636d+00,0.372006d+00, &
133 0.377447d+00,0.382966d+00,0.388567d+00,0.394233d+00,0.399981d+00, &
134 0.405806d+00,0.411714d+00,0.417699d+00,0.423772d+00,0.429914d+00, &
135 0.436145d+00,0.442468d+00,0.448862d+00,0.455359d+00,0.461930d+00, &
136 0.468596d+00,0.475348d+00,0.482186d+00,0.489124d+00,0.496160d+00, &
137 0.503278d+00,0.510497d+00,0.517808d+00,0.525224d+00,0.532737d+00, &
138 0.540355d+00,0.548059d+00,0.555886d+00,0.563797d+00,0.571825d+00, &
139 0.579952d+00,0.588198d+00,0.596545d+00,0.605000d+00,0.613572d+00, &
140 0.622255d+00,0.631059d+00,0.639962d+00,0.649003d+00,0.658144d+00, &
141 0.667414d+00,0.676815d+00,0.686317d+00,0.695956d+00,0.705728d+00, &
142 0.715622d+00,0.725641d+00,0.735799d+00,0.746082d+00,0.756495d+00, &
143 0.767052d+00,0.777741d+00,0.788576d+00,0.799549d+00,0.810656d+00, &
144 0.821914d+00,0.833314d+00,0.844854d+00,0.856555d+00,0.868415d+00/
145  data (es(ies),ies=286,380) / &
146 0.880404d+00,0.892575d+00,0.904877d+00,0.917350d+00,0.929974d+00, &
147 0.942771d+00,0.955724d+00,0.968837d+00,0.982127d+00,0.995600d+00, &
148 0.100921d+01,0.102304d+01,0.103700d+01,0.105116d+01,0.106549d+01, &
149 0.108002d+01,0.109471d+01,0.110962d+01,0.112469d+01,0.113995d+01, &
150 0.115542d+01,0.117107d+01,0.118693d+01,0.120298d+01,0.121923d+01, &
151 0.123569d+01,0.125234d+01,0.126923d+01,0.128631d+01,0.130362d+01, &
152 0.132114d+01,0.133887d+01,0.135683d+01,0.137500d+01,0.139342d+01, &
153 0.141205d+01,0.143091d+01,0.145000d+01,0.146933d+01,0.148892d+01, &
154 0.150874d+01,0.152881d+01,0.154912d+01,0.156970d+01,0.159049d+01, &
155 0.161159d+01,0.163293d+01,0.165452d+01,0.167640d+01,0.169852d+01, &
156 0.172091d+01,0.174359d+01,0.176653d+01,0.178977d+01,0.181332d+01, &
157 0.183709d+01,0.186119d+01,0.188559d+01,0.191028d+01,0.193524d+01, &
158 0.196054d+01,0.198616d+01,0.201208d+01,0.203829d+01,0.206485d+01, &
159 0.209170d+01,0.211885d+01,0.214637d+01,0.217424d+01,0.220242d+01, &
160 0.223092d+01,0.225979d+01,0.228899d+01,0.231855d+01,0.234845d+01, &
161 0.237874d+01,0.240937d+01,0.244040d+01,0.247176d+01,0.250349d+01, &
162 0.253560d+01,0.256814d+01,0.260099d+01,0.263431d+01,0.266800d+01, &
163 0.270207d+01,0.273656d+01,0.277145d+01,0.280671d+01,0.284248d+01, &
164 0.287859d+01,0.291516d+01,0.295219d+01,0.298962d+01,0.302746d+01/
165  data (es(ies),ies=381,475) / &
166 0.306579d+01,0.310454d+01,0.314377d+01,0.318351d+01,0.322360d+01, &
167 0.326427d+01,0.330538d+01,0.334694d+01,0.338894d+01,0.343155d+01, &
168 0.347456d+01,0.351809d+01,0.356216d+01,0.360673d+01,0.365184d+01, &
169 0.369744d+01,0.374352d+01,0.379018d+01,0.383743d+01,0.388518d+01, &
170 0.393344d+01,0.398230d+01,0.403177d+01,0.408175d+01,0.413229d+01, &
171 0.418343d+01,0.423514d+01,0.428746d+01,0.434034d+01,0.439389d+01, &
172 0.444808d+01,0.450276d+01,0.455820d+01,0.461423d+01,0.467084d+01, &
173 0.472816d+01,0.478607d+01,0.484468d+01,0.490393d+01,0.496389d+01, &
174 0.502446d+01,0.508580d+01,0.514776d+01,0.521047d+01,0.527385d+01, &
175 0.533798d+01,0.540279d+01,0.546838d+01,0.553466d+01,0.560173d+01, &
176 0.566949d+01,0.573807d+01,0.580750d+01,0.587749d+01,0.594846d+01, &
177 0.602017d+01,0.609260d+01,0.616591d+01,0.623995d+01,0.631490d+01, &
178 0.639061d+01,0.646723d+01,0.654477d+01,0.662293d+01,0.670220d+01, &
179 0.678227d+01,0.686313d+01,0.694495d+01,0.702777d+01,0.711142d+01, &
180 0.719592d+01,0.728140d+01,0.736790d+01,0.745527d+01,0.754352d+01, &
181 0.763298d+01,0.772316d+01,0.781442d+01,0.790676d+01,0.800001d+01, &
182 0.809435d+01,0.818967d+01,0.828606d+01,0.838343d+01,0.848194d+01, &
183 0.858144d+01,0.868207d+01,0.878392d+01,0.888673d+01,0.899060d+01, &
184 0.909567d+01,0.920172d+01,0.930909d+01,0.941765d+01,0.952730d+01/
185  data (es(ies),ies=476,570) / &
186 0.963821d+01,0.975022d+01,0.986352d+01,0.997793d+01,0.100937d+02, &
187 0.102105d+02,0.103287d+02,0.104481d+02,0.105688d+02,0.106909d+02, &
188 0.108143d+02,0.109387d+02,0.110647d+02,0.111921d+02,0.113207d+02, &
189 0.114508d+02,0.115821d+02,0.117149d+02,0.118490d+02,0.119847d+02, &
190 0.121216d+02,0.122601d+02,0.124002d+02,0.125416d+02,0.126846d+02, &
191 0.128290d+02,0.129747d+02,0.131224d+02,0.132712d+02,0.134220d+02, &
192 0.135742d+02,0.137278d+02,0.138831d+02,0.140403d+02,0.141989d+02, &
193 0.143589d+02,0.145211d+02,0.146845d+02,0.148501d+02,0.150172d+02, &
194 0.151858d+02,0.153564d+02,0.155288d+02,0.157029d+02,0.158786d+02, &
195 0.160562d+02,0.162358d+02,0.164174d+02,0.166004d+02,0.167858d+02, &
196 0.169728d+02,0.171620d+02,0.173528d+02,0.175455d+02,0.177406d+02, &
197 0.179372d+02,0.181363d+02,0.183372d+02,0.185400d+02,0.187453d+02, &
198 0.189523d+02,0.191613d+02,0.193728d+02,0.195866d+02,0.198024d+02, &
199 0.200200d+02,0.202401d+02,0.204626d+02,0.206871d+02,0.209140d+02, &
200 0.211430d+02,0.213744d+02,0.216085d+02,0.218446d+02,0.220828d+02, &
201 0.223241d+02,0.225671d+02,0.228132d+02,0.230615d+02,0.233120d+02, &
202 0.235651d+02,0.238211d+02,0.240794d+02,0.243404d+02,0.246042d+02, &
203 0.248704d+02,0.251390d+02,0.254109d+02,0.256847d+02,0.259620d+02, &
204 0.262418d+02,0.265240d+02,0.268092d+02,0.270975d+02,0.273883d+02/
205  data (es(ies),ies=571,665) / &
206 0.276822d+02,0.279792d+02,0.282789d+02,0.285812d+02,0.288867d+02, &
207 0.291954d+02,0.295075d+02,0.298222d+02,0.301398d+02,0.304606d+02, &
208 0.307848d+02,0.311119d+02,0.314424d+02,0.317763d+02,0.321133d+02, &
209 0.324536d+02,0.327971d+02,0.331440d+02,0.334940d+02,0.338475d+02, &
210 0.342050d+02,0.345654d+02,0.349295d+02,0.352975d+02,0.356687d+02, &
211 0.360430d+02,0.364221d+02,0.368042d+02,0.371896d+02,0.375790d+02, &
212 0.379725d+02,0.383692d+02,0.387702d+02,0.391744d+02,0.395839d+02, &
213 0.399958d+02,0.404118d+02,0.408325d+02,0.412574d+02,0.416858d+02, &
214 0.421188d+02,0.425551d+02,0.429962d+02,0.434407d+02,0.438910d+02, &
215 0.443439d+02,0.448024d+02,0.452648d+02,0.457308d+02,0.462018d+02, &
216 0.466775d+02,0.471582d+02,0.476428d+02,0.481313d+02,0.486249d+02, &
217 0.491235d+02,0.496272d+02,0.501349d+02,0.506479d+02,0.511652d+02, &
218 0.516876d+02,0.522142d+02,0.527474d+02,0.532836d+02,0.538266d+02, &
219 0.543737d+02,0.549254d+02,0.554839d+02,0.560456d+02,0.566142d+02, &
220 0.571872d+02,0.577662d+02,0.583498d+02,0.589392d+02,0.595347d+02, &
221 0.601346d+02,0.607410d+02,0.613519d+02,0.619689d+02,0.625922d+02, &
222 0.632204d+02,0.638550d+02,0.644959d+02,0.651418d+02,0.657942d+02, &
223 0.664516d+02,0.671158d+02,0.677864d+02,0.684624d+02,0.691451d+02, &
224 0.698345d+02,0.705293d+02,0.712312d+02,0.719398d+02,0.726542d+02/
225  data (es(ies),ies=666,760) / &
226 0.733754d+02,0.741022d+02,0.748363d+02,0.755777d+02,0.763247d+02, &
227 0.770791d+02,0.778394d+02,0.786088d+02,0.793824d+02,0.801653d+02, &
228 0.809542d+02,0.817509d+02,0.825536d+02,0.833643d+02,0.841828d+02, &
229 0.850076d+02,0.858405d+02,0.866797d+02,0.875289d+02,0.883827d+02, &
230 0.892467d+02,0.901172d+02,0.909962d+02,0.918818d+02,0.927760d+02, &
231 0.936790d+02,0.945887d+02,0.955071d+02,0.964346d+02,0.973689d+02, &
232 0.983123d+02,0.992648d+02,0.100224d+03,0.101193d+03,0.102169d+03, &
233 0.103155d+03,0.104150d+03,0.105152d+03,0.106164d+03,0.107186d+03, &
234 0.108217d+03,0.109256d+03,0.110303d+03,0.111362d+03,0.112429d+03, &
235 0.113503d+03,0.114588d+03,0.115684d+03,0.116789d+03,0.117903d+03, &
236 0.119028d+03,0.120160d+03,0.121306d+03,0.122460d+03,0.123623d+03, &
237 0.124796d+03,0.125981d+03,0.127174d+03,0.128381d+03,0.129594d+03, &
238 0.130822d+03,0.132058d+03,0.133306d+03,0.134563d+03,0.135828d+03, &
239 0.137109d+03,0.138402d+03,0.139700d+03,0.141017d+03,0.142338d+03, &
240 0.143676d+03,0.145025d+03,0.146382d+03,0.147753d+03,0.149133d+03, &
241 0.150529d+03,0.151935d+03,0.153351d+03,0.154783d+03,0.156222d+03, &
242 0.157678d+03,0.159148d+03,0.160624d+03,0.162117d+03,0.163621d+03, &
243 0.165142d+03,0.166674d+03,0.168212d+03,0.169772d+03,0.171340d+03, &
244 0.172921d+03,0.174522d+03,0.176129d+03,0.177755d+03,0.179388d+03/
245  data (es(ies),ies=761,855) / &
246 0.181040d+03,0.182707d+03,0.184382d+03,0.186076d+03,0.187782d+03, &
247 0.189503d+03,0.191240d+03,0.192989d+03,0.194758d+03,0.196535d+03, &
248 0.198332d+03,0.200141d+03,0.201963d+03,0.203805d+03,0.205656d+03, &
249 0.207532d+03,0.209416d+03,0.211317d+03,0.213236d+03,0.215167d+03, &
250 0.217121d+03,0.219087d+03,0.221067d+03,0.223064d+03,0.225080d+03, &
251 0.227113d+03,0.229160d+03,0.231221d+03,0.233305d+03,0.235403d+03, &
252 0.237520d+03,0.239655d+03,0.241805d+03,0.243979d+03,0.246163d+03, &
253 0.248365d+03,0.250593d+03,0.252830d+03,0.255093d+03,0.257364d+03, &
254 0.259667d+03,0.261979d+03,0.264312d+03,0.266666d+03,0.269034d+03, &
255 0.271430d+03,0.273841d+03,0.276268d+03,0.278722d+03,0.281185d+03, &
256 0.283677d+03,0.286190d+03,0.288714d+03,0.291266d+03,0.293834d+03, &
257 0.296431d+03,0.299045d+03,0.301676d+03,0.304329d+03,0.307006d+03, &
258 0.309706d+03,0.312423d+03,0.315165d+03,0.317930d+03,0.320705d+03, &
259 0.323519d+03,0.326350d+03,0.329199d+03,0.332073d+03,0.334973d+03, &
260 0.337897d+03,0.340839d+03,0.343800d+03,0.346794d+03,0.349806d+03, &
261 0.352845d+03,0.355918d+03,0.358994d+03,0.362112d+03,0.365242d+03, &
262 0.368407d+03,0.371599d+03,0.374802d+03,0.378042d+03,0.381293d+03, &
263 0.384588d+03,0.387904d+03,0.391239d+03,0.394604d+03,0.397988d+03, &
264 0.401411d+03,0.404862d+03,0.408326d+03,0.411829d+03,0.415352d+03/
265  data (es(ies),ies=856,950) / &
266 0.418906d+03,0.422490d+03,0.426095d+03,0.429740d+03,0.433398d+03, &
267 0.437097d+03,0.440827d+03,0.444570d+03,0.448354d+03,0.452160d+03, &
268 0.455999d+03,0.459870d+03,0.463765d+03,0.467702d+03,0.471652d+03, &
269 0.475646d+03,0.479674d+03,0.483715d+03,0.487811d+03,0.491911d+03, &
270 0.496065d+03,0.500244d+03,0.504448d+03,0.508698d+03,0.512961d+03, &
271 0.517282d+03,0.521617d+03,0.525989d+03,0.530397d+03,0.534831d+03, &
272 0.539313d+03,0.543821d+03,0.548355d+03,0.552938d+03,0.557549d+03, &
273 0.562197d+03,0.566884d+03,0.571598d+03,0.576351d+03,0.581131d+03, &
274 0.585963d+03,0.590835d+03,0.595722d+03,0.600663d+03,0.605631d+03, &
275 0.610641d+03,0.615151d+03,0.619625d+03,0.624140d+03,0.628671d+03, &
276 0.633243d+03,0.637845d+03,0.642465d+03,0.647126d+03,0.651806d+03, &
277 0.656527d+03,0.661279d+03,0.666049d+03,0.670861d+03,0.675692d+03, &
278 0.680566d+03,0.685471d+03,0.690396d+03,0.695363d+03,0.700350d+03, &
279 0.705381d+03,0.710444d+03,0.715527d+03,0.720654d+03,0.725801d+03, &
280 0.730994d+03,0.736219d+03,0.741465d+03,0.746756d+03,0.752068d+03, &
281 0.757426d+03,0.762819d+03,0.768231d+03,0.773692d+03,0.779172d+03, &
282 0.784701d+03,0.790265d+03,0.795849d+03,0.801483d+03,0.807137d+03, &
283 0.812842d+03,0.818582d+03,0.824343d+03,0.830153d+03,0.835987d+03, &
284 0.841871d+03,0.847791d+03,0.853733d+03,0.859727d+03,0.865743d+03/
285  data (es(ies),ies=951,1045) / &
286 0.871812d+03,0.877918d+03,0.884046d+03,0.890228d+03,0.896433d+03, &
287 0.902690d+03,0.908987d+03,0.915307d+03,0.921681d+03,0.928078d+03, &
288 0.934531d+03,0.941023d+03,0.947539d+03,0.954112d+03,0.960708d+03, &
289 0.967361d+03,0.974053d+03,0.980771d+03,0.987545d+03,0.994345d+03, &
290 0.100120d+04,0.100810d+04,0.101502d+04,0.102201d+04,0.102902d+04, &
291 0.103608d+04,0.104320d+04,0.105033d+04,0.105753d+04,0.106475d+04, &
292 0.107204d+04,0.107936d+04,0.108672d+04,0.109414d+04,0.110158d+04, &
293 0.110908d+04,0.111663d+04,0.112421d+04,0.113185d+04,0.113952d+04, &
294 0.114725d+04,0.115503d+04,0.116284d+04,0.117071d+04,0.117861d+04, &
295 0.118658d+04,0.119459d+04,0.120264d+04,0.121074d+04,0.121888d+04, &
296 0.122709d+04,0.123534d+04,0.124362d+04,0.125198d+04,0.126036d+04, &
297 0.126881d+04,0.127731d+04,0.128584d+04,0.129444d+04,0.130307d+04, &
298 0.131177d+04,0.132053d+04,0.132931d+04,0.133817d+04,0.134705d+04, &
299 0.135602d+04,0.136503d+04,0.137407d+04,0.138319d+04,0.139234d+04, &
300 0.140156d+04,0.141084d+04,0.142015d+04,0.142954d+04,0.143896d+04, &
301 0.144845d+04,0.145800d+04,0.146759d+04,0.147725d+04,0.148694d+04, &
302 0.149672d+04,0.150655d+04,0.151641d+04,0.152635d+04,0.153633d+04, &
303 0.154639d+04,0.155650d+04,0.156665d+04,0.157688d+04,0.158715d+04, &
304 0.159750d+04,0.160791d+04,0.161836d+04,0.162888d+04,0.163945d+04/
305  data (es(ies),ies=1046,1140) / &
306 0.165010d+04,0.166081d+04,0.167155d+04,0.168238d+04,0.169325d+04, &
307 0.170420d+04,0.171522d+04,0.172627d+04,0.173741d+04,0.174859d+04, &
308 0.175986d+04,0.177119d+04,0.178256d+04,0.179402d+04,0.180552d+04, &
309 0.181711d+04,0.182877d+04,0.184046d+04,0.185224d+04,0.186407d+04, &
310 0.187599d+04,0.188797d+04,0.190000d+04,0.191212d+04,0.192428d+04, &
311 0.193653d+04,0.194886d+04,0.196122d+04,0.197368d+04,0.198618d+04, &
312 0.199878d+04,0.201145d+04,0.202416d+04,0.203698d+04,0.204983d+04, &
313 0.206278d+04,0.207580d+04,0.208887d+04,0.210204d+04,0.211525d+04, &
314 0.212856d+04,0.214195d+04,0.215538d+04,0.216892d+04,0.218249d+04, &
315 0.219618d+04,0.220994d+04,0.222375d+04,0.223766d+04,0.225161d+04, &
316 0.226567d+04,0.227981d+04,0.229399d+04,0.230829d+04,0.232263d+04, &
317 0.233708d+04,0.235161d+04,0.236618d+04,0.238087d+04,0.239560d+04, &
318 0.241044d+04,0.242538d+04,0.244035d+04,0.245544d+04,0.247057d+04, &
319 0.248583d+04,0.250116d+04,0.251654d+04,0.253204d+04,0.254759d+04, &
320 0.256325d+04,0.257901d+04,0.259480d+04,0.261073d+04,0.262670d+04, &
321 0.264279d+04,0.265896d+04,0.267519d+04,0.269154d+04,0.270794d+04, &
322 0.272447d+04,0.274108d+04,0.275774d+04,0.277453d+04,0.279137d+04, &
323 0.280834d+04,0.282540d+04,0.284251d+04,0.285975d+04,0.287704d+04, &
324 0.289446d+04,0.291198d+04,0.292954d+04,0.294725d+04,0.296499d+04/
325  data (es(ies),ies=1141,1235) / &
326 0.298288d+04,0.300087d+04,0.301890d+04,0.303707d+04,0.305529d+04, &
327 0.307365d+04,0.309211d+04,0.311062d+04,0.312927d+04,0.314798d+04, &
328 0.316682d+04,0.318577d+04,0.320477d+04,0.322391d+04,0.324310d+04, &
329 0.326245d+04,0.328189d+04,0.330138d+04,0.332103d+04,0.334073d+04, &
330 0.336058d+04,0.338053d+04,0.340054d+04,0.342069d+04,0.344090d+04, &
331 0.346127d+04,0.348174d+04,0.350227d+04,0.352295d+04,0.354369d+04, &
332 0.356458d+04,0.358559d+04,0.360664d+04,0.362787d+04,0.364914d+04, &
333 0.367058d+04,0.369212d+04,0.371373d+04,0.373548d+04,0.375731d+04, &
334 0.377929d+04,0.380139d+04,0.382355d+04,0.384588d+04,0.386826d+04, &
335 0.389081d+04,0.391348d+04,0.393620d+04,0.395910d+04,0.398205d+04, &
336 0.400518d+04,0.402843d+04,0.405173d+04,0.407520d+04,0.409875d+04, &
337 0.412246d+04,0.414630d+04,0.417019d+04,0.419427d+04,0.421840d+04, &
338 0.424272d+04,0.426715d+04,0.429165d+04,0.431634d+04,0.434108d+04, &
339 0.436602d+04,0.439107d+04,0.441618d+04,0.444149d+04,0.446685d+04, &
340 0.449241d+04,0.451810d+04,0.454385d+04,0.456977d+04,0.459578d+04, &
341 0.462197d+04,0.464830d+04,0.467468d+04,0.470127d+04,0.472792d+04, &
342 0.475477d+04,0.478175d+04,0.480880d+04,0.483605d+04,0.486336d+04, &
343 0.489087d+04,0.491853d+04,0.494623d+04,0.497415d+04,0.500215d+04, &
344 0.503034d+04,0.505867d+04,0.508707d+04,0.511568d+04,0.514436d+04/
345  data (es(ies),ies=1236,1330) / &
346 0.517325d+04,0.520227d+04,0.523137d+04,0.526068d+04,0.529005d+04, &
347 0.531965d+04,0.534939d+04,0.537921d+04,0.540923d+04,0.543932d+04, &
348 0.546965d+04,0.550011d+04,0.553064d+04,0.556139d+04,0.559223d+04, &
349 0.562329d+04,0.565449d+04,0.568577d+04,0.571727d+04,0.574884d+04, &
350 0.578064d+04,0.581261d+04,0.584464d+04,0.587692d+04,0.590924d+04, &
351 0.594182d+04,0.597455d+04,0.600736d+04,0.604039d+04,0.607350d+04, &
352 0.610685d+04,0.614036d+04,0.617394d+04,0.620777d+04,0.624169d+04, &
353 0.627584d+04,0.631014d+04,0.634454d+04,0.637918d+04,0.641390d+04, &
354 0.644887d+04,0.648400d+04,0.651919d+04,0.655467d+04,0.659021d+04, &
355 0.662599d+04,0.666197d+04,0.669800d+04,0.673429d+04,0.677069d+04, &
356 0.680735d+04,0.684415d+04,0.688104d+04,0.691819d+04,0.695543d+04, &
357 0.699292d+04,0.703061d+04,0.706837d+04,0.710639d+04,0.714451d+04, &
358 0.718289d+04,0.722143d+04,0.726009d+04,0.729903d+04,0.733802d+04, &
359 0.737729d+04,0.741676d+04,0.745631d+04,0.749612d+04,0.753602d+04, &
360 0.757622d+04,0.761659d+04,0.765705d+04,0.769780d+04,0.773863d+04, &
361 0.777975d+04,0.782106d+04,0.786246d+04,0.790412d+04,0.794593d+04, &
362 0.798802d+04,0.803028d+04,0.807259d+04,0.811525d+04,0.815798d+04, &
363 0.820102d+04,0.824427d+04,0.828757d+04,0.833120d+04,0.837493d+04, &
364 0.841895d+04,0.846313d+04,0.850744d+04,0.855208d+04,0.859678d+04/
365  data (es(ies),ies=1331,1425) / &
366 0.864179d+04,0.868705d+04,0.873237d+04,0.877800d+04,0.882374d+04, &
367 0.886979d+04,0.891603d+04,0.896237d+04,0.900904d+04,0.905579d+04, &
368 0.910288d+04,0.915018d+04,0.919758d+04,0.924529d+04,0.929310d+04, &
369 0.934122d+04,0.938959d+04,0.943804d+04,0.948687d+04,0.953575d+04, &
370 0.958494d+04,0.963442d+04,0.968395d+04,0.973384d+04,0.978383d+04, &
371 0.983412d+04,0.988468d+04,0.993534d+04,0.998630d+04,0.100374d+05, &
372 0.100888d+05,0.101406d+05,0.101923d+05,0.102444d+05,0.102966d+05, &
373 0.103492d+05,0.104020d+05,0.104550d+05,0.105082d+05,0.105616d+05, &
374 0.106153d+05,0.106693d+05,0.107234d+05,0.107779d+05,0.108325d+05, &
375 0.108874d+05,0.109425d+05,0.109978d+05,0.110535d+05,0.111092d+05, &
376 0.111653d+05,0.112217d+05,0.112782d+05,0.113350d+05,0.113920d+05, &
377 0.114493d+05,0.115070d+05,0.115646d+05,0.116228d+05,0.116809d+05, &
378 0.117396d+05,0.117984d+05,0.118574d+05,0.119167d+05,0.119762d+05, &
379 0.120360d+05,0.120962d+05,0.121564d+05,0.122170d+05,0.122778d+05, &
380 0.123389d+05,0.124004d+05,0.124619d+05,0.125238d+05,0.125859d+05, &
381 0.126484d+05,0.127111d+05,0.127739d+05,0.128372d+05,0.129006d+05, &
382 0.129644d+05,0.130285d+05,0.130927d+05,0.131573d+05,0.132220d+05, &
383 0.132872d+05,0.133526d+05,0.134182d+05,0.134842d+05,0.135503d+05, &
384 0.136168d+05,0.136836d+05,0.137505d+05,0.138180d+05,0.138854d+05/
385  data (es(ies),ies=1426,1520) / &
386 0.139534d+05,0.140216d+05,0.140900d+05,0.141588d+05,0.142277d+05, &
387 0.142971d+05,0.143668d+05,0.144366d+05,0.145069d+05,0.145773d+05, &
388 0.146481d+05,0.147192d+05,0.147905d+05,0.148622d+05,0.149341d+05, &
389 0.150064d+05,0.150790d+05,0.151517d+05,0.152250d+05,0.152983d+05, &
390 0.153721d+05,0.154462d+05,0.155205d+05,0.155952d+05,0.156701d+05, &
391 0.157454d+05,0.158211d+05,0.158969d+05,0.159732d+05,0.160496d+05, &
392 0.161265d+05,0.162037d+05,0.162811d+05,0.163589d+05,0.164369d+05, &
393 0.165154d+05,0.165942d+05,0.166732d+05,0.167526d+05,0.168322d+05, &
394 0.169123d+05,0.169927d+05,0.170733d+05,0.171543d+05,0.172356d+05, &
395 0.173173d+05,0.173993d+05,0.174815d+05,0.175643d+05,0.176471d+05, &
396 0.177305d+05,0.178143d+05,0.178981d+05,0.179826d+05,0.180671d+05, &
397 0.181522d+05,0.182377d+05,0.183232d+05,0.184093d+05,0.184955d+05, &
398 0.185823d+05,0.186695d+05,0.187568d+05,0.188447d+05,0.189326d+05, &
399 0.190212d+05,0.191101d+05,0.191991d+05,0.192887d+05,0.193785d+05, &
400 0.194688d+05,0.195595d+05,0.196503d+05,0.197417d+05,0.198332d+05, &
401 0.199253d+05,0.200178d+05,0.201105d+05,0.202036d+05,0.202971d+05, &
402 0.203910d+05,0.204853d+05,0.205798d+05,0.206749d+05,0.207701d+05, &
403 0.208659d+05,0.209621d+05,0.210584d+05,0.211554d+05,0.212524d+05, &
404 0.213501d+05,0.214482d+05,0.215465d+05,0.216452d+05,0.217442d+05/
405  data (es(ies),ies=1521,1552) / &
406 0.218439d+05,0.219439d+05,0.220440d+05,0.221449d+05,0.222457d+05, &
407 0.223473d+05,0.224494d+05,0.225514d+05,0.226542d+05,0.227571d+05, &
408 0.228606d+05,0.229646d+05,0.230687d+05,0.231734d+05,0.232783d+05, &
409 0.233839d+05,0.234898d+05,0.235960d+05,0.237027d+05,0.238097d+05, &
410 0.239173d+05,0.240254d+05,0.241335d+05,0.242424d+05,0.243514d+05, &
411 0.244611d+05,0.245712d+05,0.246814d+05,0.247923d+05,0.249034d+05, &
412 0.250152d+05,0.250152d+05/
413 
414 do i = 1, npnts
415 
416  ! Compute the factor that converts from sat vapour pressure in a
417  ! pure water system to sat vapour pressure in air, FSUBW. This formula
418  ! is taken from equation A4.7 of Adrian Gill's book: Atmosphere-Ocean
419  ! Dynamics. Note that his formula works in terms of pressure in MB and
420  ! temperature in Celsius, so conversion of units leads to the slightly
421  ! different equation used here.
422 
423  fsubw = one + 1.0e-8_kind_real * p(i) * &
424  (4.5_kind_real+ 6.0e-4_kind_real * (t(i) - t0c) * (t(i) - t0c))
425 
426  ! use the lookup table to find saturated vapour pressure, and stored it in QS.
427 
428  tt = max(t_low, t(i))
429  tt = min(t_high, tt)
430 
431  atable = (tt - t_low + delta_t) / delta_t
432  itable = atable
433  atable = atable - itable
434 
435  qs(i) = (one - atable) * es(itable) + atable * es(itable + 1)
436 
437  ! Multiply by FSUBW to convert to saturated vapour pressure in air (equation A4.6
438  ! of Adrian Gill's book)
439 
440  qs(i) = qs(i) * fsubw
441 
442  ! Now form the accurate expression for QS, which is a rearranged
443  ! version of equation A4.3 of Gill's book.
444 
445  ! Note that at very low pressures we apply a fix, to prevent a
446  ! singularity (Qsat tends to 1.0 kg/kg).
447 
448  qs(i) = (epsilon * qs(i)) / (max(p(i), qs(i)) - one_minus_epsilon * qs(i))
449 
450 end do
451 
452 end subroutine ops_qsat
453 
454 !-------------------------------------------------------------------------------
455 !> Saturation Specific Humidity Scheme: Vapour to Liquid.
456 !!
457 !! \details Heritage: Ops_QsatWat.inc
458 !!
459 !! Returns a saturation mixing ratio given a temperature and pressure
460 !! using saturation vapour pressure calculated using the Goff-Gratch
461 !! Formulaie, adopted by the WMO as taken from Landolt-Bornstein, 1987
462 !! Numerical Data and Functional Relationships in Science and
463 !! Technology. Group V/Vol 4B Meteorology. Physical and Chemical
464 !! Properties of Air, P35
465 !!
466 !! Values in the lookup table are over water above and below 0 degrees C.
467 !!
468 !! Note: For vapour pressure oever water this formula is valid for
469 !! temperatures between 373K and 223K. The values for saturated vapour
470 !! over water in the lookup table below are out of the lower end of
471 !! this range. However it is standard WMO practice to use the formula
472 !! below its accepted range for use with the calculation of dew points
473 !! in the upper atmosphere
474 !!
475 !! Method: <br>
476 !! uses lookup tables to find eSAT, calculates qSAT directly from that.
477 !!
478 !! \author Met Office
479 !!
480 !! \date 09/06/2020: Created
481 !!
482 subroutine ops_qsatwat (QS, &
483  T, &
484  P, &
485  npnts)
486 
487 implicit none
488 
489 ! subroutine arguments:
490 integer :: npnts ! Points (=horizontal dimensions) being processed by qSAT scheme.
491 real(kind=kind_real) :: t(npnts) ! Temperature (K).
492 real(kind=kind_real) :: p(npnts) ! Pressure (Pa).
493 real(kind=kind_real) :: qs(npnts) ! Saturation mixing ratio at temperature T and pressure P (KG/KG)
494 
495 ! Local declarations:
496 real(kind=kind_real), parameter :: one_minus_epsilon = one - epsilon
497 real(kind=kind_real), parameter :: t_low = 183.15_kind_real
498 real(kind=kind_real), parameter :: t_high = 338.15_kind_real
499 real(kind=kind_real), parameter :: delta_t = 0.1_kind_real
500 integer, parameter :: n = ((t_high - t_low + (delta_t * half)) / delta_t) + one ! gives N = 1551
501 integer :: itable
502 real(kind=kind_real) :: atable
503 real(kind=kind_real) :: fsubw
504 real(kind=kind_real) :: tt
505 integer :: i
506 integer :: ies
507 real(kind=kind_real) :: es(0:n + 1) ! Table of saturation water vapour pressure (PA)
508 character(len=*), parameter :: routinename = "Ops_QsatWat"
509 
510 ! Note: 0 element is a repeat of 1st element to cater for special case
511 ! of low temperatures (.LE.T_low) for which the array index is
512 ! rounded down due to machine precision.
513 
514 data (es(ies),ies= 0, 95) / 0.186905d-01, &
515 0.186905d-01,0.190449d-01,0.194059d-01,0.197727d-01,0.201462d-01, &
516 0.205261d-01,0.209122d-01,0.213052d-01,0.217050d-01,0.221116d-01, &
517 0.225252d-01,0.229463d-01,0.233740d-01,0.238090d-01,0.242518d-01, &
518 0.247017d-01,0.251595d-01,0.256252d-01,0.260981d-01,0.265795d-01, &
519 0.270691d-01,0.275667d-01,0.280733d-01,0.285876d-01,0.291105d-01, &
520 0.296429d-01,0.301835d-01,0.307336d-01,0.312927d-01,0.318611d-01, &
521 0.324390d-01,0.330262d-01,0.336232d-01,0.342306d-01,0.348472d-01, &
522 0.354748d-01,0.361117d-01,0.367599d-01,0.374185d-01,0.380879d-01, &
523 0.387689d-01,0.394602d-01,0.401626d-01,0.408771d-01,0.416033d-01, &
524 0.423411d-01,0.430908d-01,0.438524d-01,0.446263d-01,0.454124d-01, &
525 0.462122d-01,0.470247d-01,0.478491d-01,0.486874d-01,0.495393d-01, &
526 0.504057d-01,0.512847d-01,0.521784d-01,0.530853d-01,0.540076d-01, &
527 0.549444d-01,0.558959d-01,0.568633d-01,0.578448d-01,0.588428d-01, &
528 0.598566d-01,0.608858d-01,0.619313d-01,0.629926d-01,0.640706d-01, &
529 0.651665d-01,0.662795d-01,0.674095d-01,0.685570d-01,0.697219d-01, &
530 0.709063d-01,0.721076d-01,0.733284d-01,0.745679d-01,0.758265d-01, &
531 0.771039d-01,0.784026d-01,0.797212d-01,0.810577d-01,0.824164d-01, &
532 0.837971d-01,0.851970d-01,0.866198d-01,0.880620d-01,0.895281d-01, &
533 0.910178d-01,0.925278d-01,0.940622d-01,0.956177d-01,0.971984d-01/
534 data (es(ies),ies= 96,190) / &
535 0.988051d-01,0.100433d+00,0.102085d+00,0.103764d+00,0.105467d+00, &
536 0.107196d+00,0.108953d+00,0.110732d+00,0.112541d+00,0.114376d+00, &
537 0.116238d+00,0.118130d+00,0.120046d+00,0.121993d+00,0.123969d+00, &
538 0.125973d+00,0.128009d+00,0.130075d+00,0.132167d+00,0.134296d+00, &
539 0.136452d+00,0.138642d+00,0.140861d+00,0.143115d+00,0.145404d+00, &
540 0.147723d+00,0.150078d+00,0.152466d+00,0.154889d+00,0.157346d+00, &
541 0.159841d+00,0.162372d+00,0.164939d+00,0.167545d+00,0.170185d+00, &
542 0.172866d+00,0.175584d+00,0.178340d+00,0.181139d+00,0.183977d+00, &
543 0.186855d+00,0.189773d+00,0.192737d+00,0.195736d+00,0.198783d+00, &
544 0.201875d+00,0.205007d+00,0.208186d+00,0.211409d+00,0.214676d+00, &
545 0.217993d+00,0.221355d+00,0.224764d+00,0.228220d+00,0.231728d+00, &
546 0.235284d+00,0.238888d+00,0.242542d+00,0.246251d+00,0.250010d+00, &
547 0.253821d+00,0.257688d+00,0.261602d+00,0.265575d+00,0.269607d+00, &
548 0.273689d+00,0.277830d+00,0.282027d+00,0.286287d+00,0.290598d+00, &
549 0.294972d+00,0.299405d+00,0.303904d+00,0.308462d+00,0.313082d+00, &
550 0.317763d+00,0.322512d+00,0.327324d+00,0.332201d+00,0.337141d+00, &
551 0.342154d+00,0.347234d+00,0.352387d+00,0.357601d+00,0.362889d+00, &
552 0.368257d+00,0.373685d+00,0.379194d+00,0.384773d+00,0.390433d+00, &
553 0.396159d+00,0.401968d+00,0.407861d+00,0.413820d+00,0.419866d+00/
554 data (es(ies),ies=191,285) / &
555 0.425999d+00,0.432203d+00,0.438494d+00,0.444867d+00,0.451332d+00, &
556 0.457879d+00,0.464510d+00,0.471226d+00,0.478037d+00,0.484935d+00, &
557 0.491920d+00,0.499005d+00,0.506181d+00,0.513447d+00,0.520816d+00, &
558 0.528279d+00,0.535835d+00,0.543497d+00,0.551256d+00,0.559113d+00, &
559 0.567081d+00,0.575147d+00,0.583315d+00,0.591585d+00,0.599970d+00, &
560 0.608472d+00,0.617069d+00,0.625785d+00,0.634609d+00,0.643556d+00, &
561 0.652611d+00,0.661782d+00,0.671077d+00,0.680487d+00,0.690015d+00, &
562 0.699656d+00,0.709433d+00,0.719344d+00,0.729363d+00,0.739518d+00, &
563 0.749795d+00,0.760217d+00,0.770763d+00,0.781454d+00,0.792258d+00, &
564 0.803208d+00,0.814309d+00,0.825528d+00,0.836914d+00,0.848422d+00, &
565 0.860086d+00,0.871891d+00,0.883837d+00,0.895944d+00,0.908214d+00, &
566 0.920611d+00,0.933175d+00,0.945890d+00,0.958776d+00,0.971812d+00, &
567 0.985027d+00,0.998379d+00,0.101193d+01,0.102561d+01,0.103949d+01, &
568 0.105352d+01,0.106774d+01,0.108213d+01,0.109669d+01,0.111144d+01, &
569 0.112636d+01,0.114148d+01,0.115676d+01,0.117226d+01,0.118791d+01, &
570 0.120377d+01,0.121984d+01,0.123608d+01,0.125252d+01,0.126919d+01, &
571 0.128604d+01,0.130309d+01,0.132036d+01,0.133782d+01,0.135549d+01, &
572 0.137339d+01,0.139150d+01,0.140984d+01,0.142839d+01,0.144715d+01, &
573 0.146616d+01,0.148538d+01,0.150482d+01,0.152450d+01,0.154445d+01/
574 data (es(ies),ies=286,380) / &
575 0.156459d+01,0.158502d+01,0.160564d+01,0.162654d+01,0.164766d+01, &
576 0.166906d+01,0.169070d+01,0.171257d+01,0.173473d+01,0.175718d+01, &
577 0.177984d+01,0.180282d+01,0.182602d+01,0.184951d+01,0.187327d+01, &
578 0.189733d+01,0.192165d+01,0.194629d+01,0.197118d+01,0.199636d+01, &
579 0.202185d+01,0.204762d+01,0.207372d+01,0.210010d+01,0.212678d+01, &
580 0.215379d+01,0.218109d+01,0.220873d+01,0.223668d+01,0.226497d+01, &
581 0.229357d+01,0.232249d+01,0.235176d+01,0.238134d+01,0.241129d+01, &
582 0.244157d+01,0.247217d+01,0.250316d+01,0.253447d+01,0.256617d+01, &
583 0.259821d+01,0.263064d+01,0.266341d+01,0.269661d+01,0.273009d+01, &
584 0.276403d+01,0.279834d+01,0.283302d+01,0.286811d+01,0.290358d+01, &
585 0.293943d+01,0.297571d+01,0.301236d+01,0.304946d+01,0.308702d+01, &
586 0.312491d+01,0.316326d+01,0.320208d+01,0.324130d+01,0.328092d+01, &
587 0.332102d+01,0.336162d+01,0.340264d+01,0.344407d+01,0.348601d+01, &
588 0.352838d+01,0.357118d+01,0.361449d+01,0.365834d+01,0.370264d+01, &
589 0.374737d+01,0.379265d+01,0.383839d+01,0.388469d+01,0.393144d+01, &
590 0.397876d+01,0.402656d+01,0.407492d+01,0.412378d+01,0.417313d+01, &
591 0.422306d+01,0.427359d+01,0.432454d+01,0.437617d+01,0.442834d+01, &
592 0.448102d+01,0.453433d+01,0.458816d+01,0.464253d+01,0.469764d+01, &
593 0.475321d+01,0.480942d+01,0.486629d+01,0.492372d+01,0.498173d+01/
594 data (es(ies),ies=381,475) / &
595 0.504041d+01,0.509967d+01,0.515962d+01,0.522029d+01,0.528142d+01, &
596 0.534337d+01,0.540595d+01,0.546912d+01,0.553292d+01,0.559757d+01, &
597 0.566273d+01,0.572864d+01,0.579532d+01,0.586266d+01,0.593075d+01, &
598 0.599952d+01,0.606895d+01,0.613918d+01,0.621021d+01,0.628191d+01, &
599 0.635433d+01,0.642755d+01,0.650162d+01,0.657639d+01,0.665188d+01, &
600 0.672823d+01,0.680532d+01,0.688329d+01,0.696198d+01,0.704157d+01, &
601 0.712206d+01,0.720319d+01,0.728534d+01,0.736829d+01,0.745204d+01, &
602 0.753671d+01,0.762218d+01,0.770860d+01,0.779588d+01,0.788408d+01, &
603 0.797314d+01,0.806318d+01,0.815408d+01,0.824599d+01,0.833874d+01, &
604 0.843254d+01,0.852721d+01,0.862293d+01,0.871954d+01,0.881724d+01, &
605 0.891579d+01,0.901547d+01,0.911624d+01,0.921778d+01,0.932061d+01, &
606 0.942438d+01,0.952910d+01,0.963497d+01,0.974181d+01,0.984982d+01, &
607 0.995887d+01,0.100690d+02,0.101804d+02,0.102926d+02,0.104063d+02, &
608 0.105210d+02,0.106367d+02,0.107536d+02,0.108719d+02,0.109912d+02, &
609 0.111116d+02,0.112333d+02,0.113563d+02,0.114804d+02,0.116056d+02, &
610 0.117325d+02,0.118602d+02,0.119892d+02,0.121197d+02,0.122513d+02, &
611 0.123844d+02,0.125186d+02,0.126543d+02,0.127912d+02,0.129295d+02, &
612 0.130691d+02,0.132101d+02,0.133527d+02,0.134965d+02,0.136415d+02, &
613 0.137882d+02,0.139361d+02,0.140855d+02,0.142366d+02,0.143889d+02/
614 data (es(ies),ies=476,570) / &
615 0.145429d+02,0.146982d+02,0.148552d+02,0.150135d+02,0.151735d+02, &
616 0.153349d+02,0.154979d+02,0.156624d+02,0.158286d+02,0.159965d+02, &
617 0.161659d+02,0.163367d+02,0.165094d+02,0.166838d+02,0.168597d+02, &
618 0.170375d+02,0.172168d+02,0.173979d+02,0.175806d+02,0.177651d+02, &
619 0.179513d+02,0.181394d+02,0.183293d+02,0.185210d+02,0.187146d+02, &
620 0.189098d+02,0.191066d+02,0.193059d+02,0.195065d+02,0.197095d+02, &
621 0.199142d+02,0.201206d+02,0.203291d+02,0.205397d+02,0.207522d+02, &
622 0.209664d+02,0.211831d+02,0.214013d+02,0.216221d+02,0.218448d+02, &
623 0.220692d+02,0.222959d+02,0.225250d+02,0.227559d+02,0.229887d+02, &
624 0.232239d+02,0.234614d+02,0.237014d+02,0.239428d+02,0.241872d+02, &
625 0.244335d+02,0.246824d+02,0.249332d+02,0.251860d+02,0.254419d+02, &
626 0.256993d+02,0.259600d+02,0.262225d+02,0.264873d+02,0.267552d+02, &
627 0.270248d+02,0.272970d+02,0.275719d+02,0.278497d+02,0.281295d+02, &
628 0.284117d+02,0.286965d+02,0.289843d+02,0.292743d+02,0.295671d+02, &
629 0.298624d+02,0.301605d+02,0.304616d+02,0.307650d+02,0.310708d+02, &
630 0.313803d+02,0.316915d+02,0.320064d+02,0.323238d+02,0.326437d+02, &
631 0.329666d+02,0.332928d+02,0.336215d+02,0.339534d+02,0.342885d+02, &
632 0.346263d+02,0.349666d+02,0.353109d+02,0.356572d+02,0.360076d+02, &
633 0.363606d+02,0.367164d+02,0.370757d+02,0.374383d+02,0.378038d+02/
634 data (es(ies),ies=571,665) / &
635 0.381727d+02,0.385453d+02,0.389206d+02,0.392989d+02,0.396807d+02, &
636 0.400663d+02,0.404555d+02,0.408478d+02,0.412428d+02,0.416417d+02, &
637 0.420445d+02,0.424502d+02,0.428600d+02,0.432733d+02,0.436900d+02, &
638 0.441106d+02,0.445343d+02,0.449620d+02,0.453930d+02,0.458280d+02, &
639 0.462672d+02,0.467096d+02,0.471561d+02,0.476070d+02,0.480610d+02, &
640 0.485186d+02,0.489813d+02,0.494474d+02,0.499170d+02,0.503909d+02, &
641 0.508693d+02,0.513511d+02,0.518376d+02,0.523277d+02,0.528232d+02, &
642 0.533213d+02,0.538240d+02,0.543315d+02,0.548437d+02,0.553596d+02, &
643 0.558802d+02,0.564046d+02,0.569340d+02,0.574672d+02,0.580061d+02, &
644 0.585481d+02,0.590963d+02,0.596482d+02,0.602041d+02,0.607649d+02, &
645 0.613311d+02,0.619025d+02,0.624779d+02,0.630574d+02,0.636422d+02, &
646 0.642324d+02,0.648280d+02,0.654278d+02,0.660332d+02,0.666426d+02, &
647 0.672577d+02,0.678771d+02,0.685034d+02,0.691328d+02,0.697694d+02, &
648 0.704103d+02,0.710556d+02,0.717081d+02,0.723639d+02,0.730269d+02, &
649 0.736945d+02,0.743681d+02,0.750463d+02,0.757309d+02,0.764214d+02, &
650 0.771167d+02,0.778182d+02,0.785246d+02,0.792373d+02,0.799564d+02, &
651 0.806804d+02,0.814109d+02,0.821479d+02,0.828898d+02,0.836384d+02, &
652 0.843922d+02,0.851525d+02,0.859198d+02,0.866920d+02,0.874712d+02, &
653 0.882574d+02,0.890486d+02,0.898470d+02,0.906525d+02,0.914634d+02/
654 data (es(ies),ies=666,760) / &
655 0.922814d+02,0.931048d+02,0.939356d+02,0.947736d+02,0.956171d+02, &
656 0.964681d+02,0.973246d+02,0.981907d+02,0.990605d+02,0.999399d+02, &
657 0.100825d+03,0.101718d+03,0.102617d+03,0.103523d+03,0.104438d+03, &
658 0.105358d+03,0.106287d+03,0.107221d+03,0.108166d+03,0.109115d+03, &
659 0.110074d+03,0.111039d+03,0.112012d+03,0.112992d+03,0.113981d+03, &
660 0.114978d+03,0.115981d+03,0.116993d+03,0.118013d+03,0.119041d+03, &
661 0.120077d+03,0.121122d+03,0.122173d+03,0.123234d+03,0.124301d+03, &
662 0.125377d+03,0.126463d+03,0.127556d+03,0.128657d+03,0.129769d+03, &
663 0.130889d+03,0.132017d+03,0.133152d+03,0.134299d+03,0.135453d+03, &
664 0.136614d+03,0.137786d+03,0.138967d+03,0.140158d+03,0.141356d+03, &
665 0.142565d+03,0.143781d+03,0.145010d+03,0.146247d+03,0.147491d+03, &
666 0.148746d+03,0.150011d+03,0.151284d+03,0.152571d+03,0.153862d+03, &
667 0.155168d+03,0.156481d+03,0.157805d+03,0.159137d+03,0.160478d+03, &
668 0.161832d+03,0.163198d+03,0.164569d+03,0.165958d+03,0.167348d+03, &
669 0.168757d+03,0.170174d+03,0.171599d+03,0.173037d+03,0.174483d+03, &
670 0.175944d+03,0.177414d+03,0.178892d+03,0.180387d+03,0.181886d+03, &
671 0.183402d+03,0.184930d+03,0.186463d+03,0.188012d+03,0.189571d+03, &
672 0.191146d+03,0.192730d+03,0.194320d+03,0.195930d+03,0.197546d+03, &
673 0.199175d+03,0.200821d+03,0.202473d+03,0.204142d+03,0.205817d+03/
674 data (es(ies),ies=761,855) / &
675 0.207510d+03,0.209216d+03,0.210928d+03,0.212658d+03,0.214398d+03, &
676 0.216152d+03,0.217920d+03,0.219698d+03,0.221495d+03,0.223297d+03, &
677 0.225119d+03,0.226951d+03,0.228793d+03,0.230654d+03,0.232522d+03, &
678 0.234413d+03,0.236311d+03,0.238223d+03,0.240151d+03,0.242090d+03, &
679 0.244049d+03,0.246019d+03,0.248000d+03,0.249996d+03,0.252009d+03, &
680 0.254037d+03,0.256077d+03,0.258128d+03,0.260200d+03,0.262284d+03, &
681 0.264384d+03,0.266500d+03,0.268629d+03,0.270779d+03,0.272936d+03, &
682 0.275110d+03,0.277306d+03,0.279509d+03,0.281734d+03,0.283966d+03, &
683 0.286227d+03,0.288494d+03,0.290780d+03,0.293083d+03,0.295398d+03, &
684 0.297737d+03,0.300089d+03,0.302453d+03,0.304841d+03,0.307237d+03, &
685 0.309656d+03,0.312095d+03,0.314541d+03,0.317012d+03,0.319496d+03, &
686 0.322005d+03,0.324527d+03,0.327063d+03,0.329618d+03,0.332193d+03, &
687 0.334788d+03,0.337396d+03,0.340025d+03,0.342673d+03,0.345329d+03, &
688 0.348019d+03,0.350722d+03,0.353440d+03,0.356178d+03,0.358938d+03, &
689 0.361718d+03,0.364513d+03,0.367322d+03,0.370160d+03,0.373012d+03, &
690 0.375885d+03,0.378788d+03,0.381691d+03,0.384631d+03,0.387579d+03, &
691 0.390556d+03,0.393556d+03,0.396563d+03,0.399601d+03,0.402646d+03, &
692 0.405730d+03,0.408829d+03,0.411944d+03,0.415083d+03,0.418236d+03, &
693 0.421422d+03,0.424632d+03,0.427849d+03,0.431099d+03,0.434365d+03/
694 data (es(ies),ies=856,950) / &
695 0.437655d+03,0.440970d+03,0.444301d+03,0.447666d+03,0.451038d+03, &
696 0.454445d+03,0.457876d+03,0.461316d+03,0.464790d+03,0.468281d+03, &
697 0.471798d+03,0.475342d+03,0.478902d+03,0.482497d+03,0.486101d+03, &
698 0.489741d+03,0.493408d+03,0.497083d+03,0.500804d+03,0.504524d+03, &
699 0.508290d+03,0.512074d+03,0.515877d+03,0.519717d+03,0.523566d+03, &
700 0.527462d+03,0.531367d+03,0.535301d+03,0.539264d+03,0.543245d+03, &
701 0.547265d+03,0.551305d+03,0.555363d+03,0.559462d+03,0.563579d+03, &
702 0.567727d+03,0.571905d+03,0.576102d+03,0.580329d+03,0.584576d+03, &
703 0.588865d+03,0.593185d+03,0.597514d+03,0.601885d+03,0.606276d+03, &
704 0.610699d+03,0.615151d+03,0.619625d+03,0.624140d+03,0.628671d+03, &
705 0.633243d+03,0.637845d+03,0.642465d+03,0.647126d+03,0.651806d+03, &
706 0.656527d+03,0.661279d+03,0.666049d+03,0.670861d+03,0.675692d+03, &
707 0.680566d+03,0.685471d+03,0.690396d+03,0.695363d+03,0.700350d+03, &
708 0.705381d+03,0.710444d+03,0.715527d+03,0.720654d+03,0.725801d+03, &
709 0.730994d+03,0.736219d+03,0.741465d+03,0.746756d+03,0.752068d+03, &
710 0.757426d+03,0.762819d+03,0.768231d+03,0.773692d+03,0.779172d+03, &
711 0.784701d+03,0.790265d+03,0.795849d+03,0.801483d+03,0.807137d+03, &
712 0.812842d+03,0.818582d+03,0.824343d+03,0.830153d+03,0.835987d+03, &
713 0.841871d+03,0.847791d+03,0.853733d+03,0.859727d+03,0.865743d+03/
714 data (es(ies),ies=951,1045) / &
715 0.871812d+03,0.877918d+03,0.884046d+03,0.890228d+03,0.896433d+03, &
716 0.902690d+03,0.908987d+03,0.915307d+03,0.921681d+03,0.928078d+03, &
717 0.934531d+03,0.941023d+03,0.947539d+03,0.954112d+03,0.960708d+03, &
718 0.967361d+03,0.974053d+03,0.980771d+03,0.987545d+03,0.994345d+03, &
719 0.100120d+04,0.100810d+04,0.101502d+04,0.102201d+04,0.102902d+04, &
720 0.103608d+04,0.104320d+04,0.105033d+04,0.105753d+04,0.106475d+04, &
721 0.107204d+04,0.107936d+04,0.108672d+04,0.109414d+04,0.110158d+04, &
722 0.110908d+04,0.111663d+04,0.112421d+04,0.113185d+04,0.113952d+04, &
723 0.114725d+04,0.115503d+04,0.116284d+04,0.117071d+04,0.117861d+04, &
724 0.118658d+04,0.119459d+04,0.120264d+04,0.121074d+04,0.121888d+04, &
725 0.122709d+04,0.123534d+04,0.124362d+04,0.125198d+04,0.126036d+04, &
726 0.126881d+04,0.127731d+04,0.128584d+04,0.129444d+04,0.130307d+04, &
727 0.131177d+04,0.132053d+04,0.132931d+04,0.133817d+04,0.134705d+04, &
728 0.135602d+04,0.136503d+04,0.137407d+04,0.138319d+04,0.139234d+04, &
729 0.140156d+04,0.141084d+04,0.142015d+04,0.142954d+04,0.143896d+04, &
730 0.144845d+04,0.145800d+04,0.146759d+04,0.147725d+04,0.148694d+04, &
731 0.149672d+04,0.150655d+04,0.151641d+04,0.152635d+04,0.153633d+04, &
732 0.154639d+04,0.155650d+04,0.156665d+04,0.157688d+04,0.158715d+04, &
733 0.159750d+04,0.160791d+04,0.161836d+04,0.162888d+04,0.163945d+04/
734 data (es(ies),ies=1046,1140) / &
735 0.165010d+04,0.166081d+04,0.167155d+04,0.168238d+04,0.169325d+04, &
736 0.170420d+04,0.171522d+04,0.172627d+04,0.173741d+04,0.174859d+04, &
737 0.175986d+04,0.177119d+04,0.178256d+04,0.179402d+04,0.180552d+04, &
738 0.181711d+04,0.182877d+04,0.184046d+04,0.185224d+04,0.186407d+04, &
739 0.187599d+04,0.188797d+04,0.190000d+04,0.191212d+04,0.192428d+04, &
740 0.193653d+04,0.194886d+04,0.196122d+04,0.197368d+04,0.198618d+04, &
741 0.199878d+04,0.201145d+04,0.202416d+04,0.203698d+04,0.204983d+04, &
742 0.206278d+04,0.207580d+04,0.208887d+04,0.210204d+04,0.211525d+04, &
743 0.212856d+04,0.214195d+04,0.215538d+04,0.216892d+04,0.218249d+04, &
744 0.219618d+04,0.220994d+04,0.222375d+04,0.223766d+04,0.225161d+04, &
745 0.226567d+04,0.227981d+04,0.229399d+04,0.230829d+04,0.232263d+04, &
746 0.233708d+04,0.235161d+04,0.236618d+04,0.238087d+04,0.239560d+04, &
747 0.241044d+04,0.242538d+04,0.244035d+04,0.245544d+04,0.247057d+04, &
748 0.248583d+04,0.250116d+04,0.251654d+04,0.253204d+04,0.254759d+04, &
749 0.256325d+04,0.257901d+04,0.259480d+04,0.261073d+04,0.262670d+04, &
750 0.264279d+04,0.265896d+04,0.267519d+04,0.269154d+04,0.270794d+04, &
751 0.272447d+04,0.274108d+04,0.275774d+04,0.277453d+04,0.279137d+04, &
752 0.280834d+04,0.282540d+04,0.284251d+04,0.285975d+04,0.287704d+04, &
753 0.289446d+04,0.291198d+04,0.292954d+04,0.294725d+04,0.296499d+04/
754 data (es(ies),ies=1141,1235) / &
755 0.298288d+04,0.300087d+04,0.301890d+04,0.303707d+04,0.305529d+04, &
756 0.307365d+04,0.309211d+04,0.311062d+04,0.312927d+04,0.314798d+04, &
757 0.316682d+04,0.318577d+04,0.320477d+04,0.322391d+04,0.324310d+04, &
758 0.326245d+04,0.328189d+04,0.330138d+04,0.332103d+04,0.334073d+04, &
759 0.336058d+04,0.338053d+04,0.340054d+04,0.342069d+04,0.344090d+04, &
760 0.346127d+04,0.348174d+04,0.350227d+04,0.352295d+04,0.354369d+04, &
761 0.356458d+04,0.358559d+04,0.360664d+04,0.362787d+04,0.364914d+04, &
762 0.367058d+04,0.369212d+04,0.371373d+04,0.373548d+04,0.375731d+04, &
763 0.377929d+04,0.380139d+04,0.382355d+04,0.384588d+04,0.386826d+04, &
764 0.389081d+04,0.391348d+04,0.393620d+04,0.395910d+04,0.398205d+04, &
765 0.400518d+04,0.402843d+04,0.405173d+04,0.407520d+04,0.409875d+04, &
766 0.412246d+04,0.414630d+04,0.417019d+04,0.419427d+04,0.421840d+04, &
767 0.424272d+04,0.426715d+04,0.429165d+04,0.431634d+04,0.434108d+04, &
768 0.436602d+04,0.439107d+04,0.441618d+04,0.444149d+04,0.446685d+04, &
769 0.449241d+04,0.451810d+04,0.454385d+04,0.456977d+04,0.459578d+04, &
770 0.462197d+04,0.464830d+04,0.467468d+04,0.470127d+04,0.472792d+04, &
771 0.475477d+04,0.478175d+04,0.480880d+04,0.483605d+04,0.486336d+04, &
772 0.489087d+04,0.491853d+04,0.494623d+04,0.497415d+04,0.500215d+04, &
773 0.503034d+04,0.505867d+04,0.508707d+04,0.511568d+04,0.514436d+04/
774 data (es(ies),ies=1236,1330) / &
775 0.517325d+04,0.520227d+04,0.523137d+04,0.526068d+04,0.529005d+04, &
776 0.531965d+04,0.534939d+04,0.537921d+04,0.540923d+04,0.543932d+04, &
777 0.546965d+04,0.550011d+04,0.553064d+04,0.556139d+04,0.559223d+04, &
778 0.562329d+04,0.565449d+04,0.568577d+04,0.571727d+04,0.574884d+04, &
779 0.578064d+04,0.581261d+04,0.584464d+04,0.587692d+04,0.590924d+04, &
780 0.594182d+04,0.597455d+04,0.600736d+04,0.604039d+04,0.607350d+04, &
781 0.610685d+04,0.614036d+04,0.617394d+04,0.620777d+04,0.624169d+04, &
782 0.627584d+04,0.631014d+04,0.634454d+04,0.637918d+04,0.641390d+04, &
783 0.644887d+04,0.648400d+04,0.651919d+04,0.655467d+04,0.659021d+04, &
784 0.662599d+04,0.666197d+04,0.669800d+04,0.673429d+04,0.677069d+04, &
785 0.680735d+04,0.684415d+04,0.688104d+04,0.691819d+04,0.695543d+04, &
786 0.699292d+04,0.703061d+04,0.706837d+04,0.710639d+04,0.714451d+04, &
787 0.718289d+04,0.722143d+04,0.726009d+04,0.729903d+04,0.733802d+04, &
788 0.737729d+04,0.741676d+04,0.745631d+04,0.749612d+04,0.753602d+04, &
789 0.757622d+04,0.761659d+04,0.765705d+04,0.769780d+04,0.773863d+04, &
790 0.777975d+04,0.782106d+04,0.786246d+04,0.790412d+04,0.794593d+04, &
791 0.798802d+04,0.803028d+04,0.807259d+04,0.811525d+04,0.815798d+04, &
792 0.820102d+04,0.824427d+04,0.828757d+04,0.833120d+04,0.837493d+04, &
793 0.841895d+04,0.846313d+04,0.850744d+04,0.855208d+04,0.859678d+04/
794 data (es(ies),ies=1331,1425) / &
795 0.864179d+04,0.868705d+04,0.873237d+04,0.877800d+04,0.882374d+04, &
796 0.886979d+04,0.891603d+04,0.896237d+04,0.900904d+04,0.905579d+04, &
797 0.910288d+04,0.915018d+04,0.919758d+04,0.924529d+04,0.929310d+04, &
798 0.934122d+04,0.938959d+04,0.943804d+04,0.948687d+04,0.953575d+04, &
799 0.958494d+04,0.963442d+04,0.968395d+04,0.973384d+04,0.978383d+04, &
800 0.983412d+04,0.988468d+04,0.993534d+04,0.998630d+04,0.100374d+05, &
801 0.100888d+05,0.101406d+05,0.101923d+05,0.102444d+05,0.102966d+05, &
802 0.103492d+05,0.104020d+05,0.104550d+05,0.105082d+05,0.105616d+05, &
803 0.106153d+05,0.106693d+05,0.107234d+05,0.107779d+05,0.108325d+05, &
804 0.108874d+05,0.109425d+05,0.109978d+05,0.110535d+05,0.111092d+05, &
805 0.111653d+05,0.112217d+05,0.112782d+05,0.113350d+05,0.113920d+05, &
806 0.114493d+05,0.115070d+05,0.115646d+05,0.116228d+05,0.116809d+05, &
807 0.117396d+05,0.117984d+05,0.118574d+05,0.119167d+05,0.119762d+05, &
808 0.120360d+05,0.120962d+05,0.121564d+05,0.122170d+05,0.122778d+05, &
809 0.123389d+05,0.124004d+05,0.124619d+05,0.125238d+05,0.125859d+05, &
810 0.126484d+05,0.127111d+05,0.127739d+05,0.128372d+05,0.129006d+05, &
811 0.129644d+05,0.130285d+05,0.130927d+05,0.131573d+05,0.132220d+05, &
812 0.132872d+05,0.133526d+05,0.134182d+05,0.134842d+05,0.135503d+05, &
813 0.136168d+05,0.136836d+05,0.137505d+05,0.138180d+05,0.138854d+05/
814 data (es(ies),ies=1426,1520) / &
815 0.139534d+05,0.140216d+05,0.140900d+05,0.141588d+05,0.142277d+05, &
816 0.142971d+05,0.143668d+05,0.144366d+05,0.145069d+05,0.145773d+05, &
817 0.146481d+05,0.147192d+05,0.147905d+05,0.148622d+05,0.149341d+05, &
818 0.150064d+05,0.150790d+05,0.151517d+05,0.152250d+05,0.152983d+05, &
819 0.153721d+05,0.154462d+05,0.155205d+05,0.155952d+05,0.156701d+05, &
820 0.157454d+05,0.158211d+05,0.158969d+05,0.159732d+05,0.160496d+05, &
821 0.161265d+05,0.162037d+05,0.162811d+05,0.163589d+05,0.164369d+05, &
822 0.165154d+05,0.165942d+05,0.166732d+05,0.167526d+05,0.168322d+05, &
823 0.169123d+05,0.169927d+05,0.170733d+05,0.171543d+05,0.172356d+05, &
824 0.173173d+05,0.173993d+05,0.174815d+05,0.175643d+05,0.176471d+05, &
825 0.177305d+05,0.178143d+05,0.178981d+05,0.179826d+05,0.180671d+05, &
826 0.181522d+05,0.182377d+05,0.183232d+05,0.184093d+05,0.184955d+05, &
827 0.185823d+05,0.186695d+05,0.187568d+05,0.188447d+05,0.189326d+05, &
828 0.190212d+05,0.191101d+05,0.191991d+05,0.192887d+05,0.193785d+05, &
829 0.194688d+05,0.195595d+05,0.196503d+05,0.197417d+05,0.198332d+05, &
830 0.199253d+05,0.200178d+05,0.201105d+05,0.202036d+05,0.202971d+05, &
831 0.203910d+05,0.204853d+05,0.205798d+05,0.206749d+05,0.207701d+05, &
832 0.208659d+05,0.209621d+05,0.210584d+05,0.211554d+05,0.212524d+05, &
833 0.213501d+05,0.214482d+05,0.215465d+05,0.216452d+05,0.217442d+05/
834 data (es(ies),ies=1521,1552) / &
835 0.218439d+05,0.219439d+05,0.220440d+05,0.221449d+05,0.222457d+05, &
836 0.223473d+05,0.224494d+05,0.225514d+05,0.226542d+05,0.227571d+05, &
837 0.228606d+05,0.229646d+05,0.230687d+05,0.231734d+05,0.232783d+05, &
838 0.233839d+05,0.234898d+05,0.235960d+05,0.237027d+05,0.238097d+05, &
839 0.239173d+05,0.240254d+05,0.241335d+05,0.242424d+05,0.243514d+05, &
840 0.244611d+05,0.245712d+05,0.246814d+05,0.247923d+05,0.249034d+05, &
841 0.250152d+05,0.250152d+05/
842 
843 do i = 1, npnts
844 
845  ! Compute the factor that converts from sat vapour pressure in a
846  ! pure water system to sat vapour pressure in air, FSUBW. This formula
847  ! is taken from equation A4.7 of Adrian Gill's Book: Atmosphere-Ocean
848  ! Dynamics. Note that his formula works in terms of pressure in MB and
849  ! temperature in Celsius, so conversion of units leads to the slightly
850  ! different equation used here.
851 
852  fsubw = one + 1.0e-8_kind_real * p(i) &
853  * (4.5_kind_real + 6.0e-4_kind_real * (t(i) - t0c) * (t(i) - t0c))
854 
855  ! use the lookup table to find saturated vapour pressure, and store it in QS
856 
857  tt = max(t_low, t(i))
858  tt = min(t_high, tt)
859 
860  atable = (tt - t_low + delta_t) / delta_t
861  itable = atable
862  atable = atable - itable
863 
864  qs(i) = (one - atable) * es(itable) + atable * es(itable + 1)
865 
866  ! Multiple by FSUBW to convert to saturated vapour pressure in air
867  ! (equation A4.6 of Adrian Gill's book)
868 
869  qs(i) = qs(i) * fsubw
870 
871  ! Now form the accurate expression for QS, which is a rearranged version of
872  ! equation A4.3 of Gill's book.
873  !
874  ! Note that at very low pressures we apply a fix, to prevent a singularity
875  ! (Qsat tends to 1.0 kg/kg)
876 
877  qs(i) = (epsilon * qs(i)) / (max(p(i), qs(i)) - one_minus_epsilon * qs(i))
878 
879 end do
880 
881 end subroutine ops_qsatwat
882 
883 !-------------------------------------------------------------------------------
884 !> Split the humidity into water vapour, liquid water and ice.
885 !!
886 !! \details Heritage: Ops_SatRad_Qsplit.f90
887 !!
888 !! if output_type=1 : Split total water content (qtotal) into <br>
889 !! water vapor content (q) and <br>
890 !! cloud liquid water content (ql) and <br>
891 !! cloud ice water content (qi) <br>
892 !!
893 !! if output_type ne 1 : Compute derivatives: (q) =dq/dqtotal <br>
894 !! (ql)=dql/dqtotal <br>
895 !! (qi)=dqi/dqtotal <br>
896 !!
897 !! \warning The derivatives are not valid if LtemperatureVar=.true. since
898 !! qsaturated depends on temperature.
899 !!
900 !! The partitioning of the excess moisture between ice and clw uses a temperature
901 !! based parametrization based on aircraft data Ref: Jones DC Reading phdthesis
902 !! p126
903 !!
904 !! \author Met Office
905 !!
906 !! \date 09/06/2020: Created
907 !!
908 subroutine ops_satrad_qsplit ( &
909  output_type, &
910  p, &
911  t, &
912  qtotal, &
913  q, &
914  ql, &
915  qi, &
916  UseQtSplitRain)
917 
918  implicit none
919 
920  ! Subroutine arguments:
921  integer, intent(in) :: output_type
922  real(kind_real), intent(in) :: p(:)
923  real(kind_real), intent(in) :: t(:)
924  real(kind_real), intent(in) :: qtotal(:)
925  real(kind_real), intent(out) :: q(size(qtotal)) ! humidity component q
926  real(kind_real), intent(out) :: ql(size(qtotal)) ! liquid component ql
927  real(kind_real), intent(out) :: qi(size(qtotal)) ! ice component qi
928  logical, intent(in) :: useqtsplitrain
929 
930  ! Local declarations:
931  integer :: i
932  real(kind_real), parameter :: lower_rh = 0.95_kind_real
933  real(kind_real), parameter :: upper_rh = 1.05_kind_real
934  real(kind_real), parameter :: split_factor = half
935  real(kind_real), parameter :: mintempql = 233.15_kind_real ! temperature (K) below which all cloud is ice
936  real(kind_real) :: qsaturated(size(qtotal))
937  real(kind_real) :: rh_qtotal(size(qtotal))
938  real(kind_real) :: qnv(size(qtotal)) ! non vapour component
939  real(kind_real) :: qc(size(qtotal)) ! cloud component
940  real(kind_real) :: v1(size(qtotal))
941  real(kind_real) :: v2(size(qtotal))
942  real(kind_real) :: v1zero
943  real(kind_real) :: v2zero
944  real(kind_real) :: w(size(qtotal))
945  real(kind_real) :: y1,y2,y3,y4
946  real(kind_real) :: intconst
947  real(kind_real) :: aconst
948  real(kind_real) :: bconst
949  real(kind_real) :: cconst
950  real(kind_real) :: dconst
951  real(kind_real) :: denom
952  real(kind_real) :: smallvalue
953  real(kind_real) :: lf(size(qtotal)) ! fraction of ql to ql+qi
954  character(len=*), parameter :: routinename = "Ops_SatRad_Qsplit"
955 
956  real(kind_real) :: qsplitrainparama !Parameters used to define proportion of
957  real(kind_real) :: qsplitrainparamb !qt that is partitioned into a rain compnenent
958  real(kind_real) :: qsplitrainparamc !
959 
960  !Qsplit_MixPhaseParam = 1 !Jones' method as default
961  qsplitrainparama = 0.15_kind_real
962  qsplitrainparamb = 0.09_kind_real
963  qsplitrainparamc = 50.0_kind_real
964 
965  ! Compute saturated water vapor profile for nlevels_q only
966 
967  call ops_qsat (qsaturated(:), & ! out
968  t(:), & ! in
969  p(:), & ! in
970  size(qtotal)) ! in
971 
972  smallvalue = one / 8.5_kind_real
973  denom = smallvalue * (upper_rh - lower_rh)
974 
975  ! don't let rh exceed 2.0 to avoid cosh function blowing up
976  rh_qtotal(:) = min(qtotal(:) / qsaturated(:), two)
977 
978  v1(:) = (rh_qtotal(:) - lower_rh) / denom
979  v2(:) = (rh_qtotal(:) - upper_rh) / denom
980 
981  y1 = one
982  y2 = split_factor
983  y3 = zero
984  y4 = split_factor
985 
986  aconst = (y2 - y1) / two
987  bconst = -(y4 - y3) / two
988  cconst = (y2 + y1) / two
989  dconst = -(y4 + y3) / two
990 
991  ! Compute fraction of ql to ql+qi based on temperature profile
992 
993  where (t(:) - zerodegc >= -0.01_kind_real) ! -0.01degc and above
994 
995  ! all ql
996  lf(:) = one
997 
998  end where
999 
1000  where (t(:) <= mintempql)
1001 
1002  ! all qi
1003  lf(:) = zero
1004 
1005  end where
1006 
1007  where (t(:) > mintempql .and. t(:) - zerodegc < -0.01_kind_real)
1008 
1009  ! Jones' parametrization
1010  lf(:) = sqrt(-1.0_kind_real * log(-0.025_kind_real * (t(:) - zerodegc)) / 70.0_kind_real)
1011 
1012  end where
1013 
1014  ! finally set LF to 0.0_kind_real for the rttov levels on which clw jacobians are not
1015  ! calculated since nlevels_mwclw < nlevels_q
1016 
1017  v1zero = -1.0_kind_real * lower_rh / denom
1018  v2zero = -1.0_kind_real * upper_rh / denom
1019  intconst = -(aconst * denom * log(cosh(v1zero)) + bconst * denom * log(cosh(v2zero)))
1020  w(:) = aconst * denom * log(cosh(v1(:))) + bconst * denom * log(cosh(v2(:))) + &
1021  (cconst + dconst) * rh_qtotal(:) + intconst
1022 
1023  ! store the components of qtotal
1024  ! ensuring that they are above lower limits
1025 
1026  if (useqtsplitrain) then
1027 
1028  ! Split qtotal into q and qnv (non-vapour part - includes
1029  ! ql, qi, qr)
1030 
1031  ! Split qtotal into q, ql ,qi
1032 
1033  do i = 1, size(qtotal)
1034 
1035  q(i) = max(w(i) * qsaturated(i), min_q)
1036  qnv(i) = max(qtotal(i) - q(i), zero)
1037 
1038  ! Split qnv into a cloud and precipitation part
1039 
1040  qc(i) = max(qsplitrainparama * (qsplitrainparamb - (qsplitrainparamb / ((qsplitrainparamc * qnv(i)) + 1))), zero)
1041 
1042  ! Finally split non-precip part into liquid and ice
1043 
1044  ql(i) = max(lf(i) * qc(i), zero)
1045  qi(i) = max((one - lf(i)) * (qc(i)), zero)
1046 
1047  end do
1048 
1049  else
1050  do i = 1, size(qtotal)
1051 
1052  q(i) = max(w(i) * qsaturated(i), min_q)
1053  ql(i) = max(lf(i) * (qtotal(i) - q(i)), zero)
1054  qi(i) = max((one - lf(i)) * (qtotal(i) - q(i)), zero)
1055 
1056  end do
1057 
1058  end if
1059 
1060  ! Values of q, ql and qi are overwritten if output_type /= 1
1061  ! and replaced with the derivatives
1062 
1063  if (output_type /= 1) then
1064 
1065  ! Compute derivates
1066  ! q = dq/dqtotal, ql = dql/dqtotal, qi=dqi/dqtotal
1067 
1068  q(:) = aconst * tanh(v1(:)) + cconst + bconst * tanh(v2(:)) + dconst
1069 
1070  if (useqtsplitrain) then
1071 
1072  ql(:) = lf(:) * qsplitrainparama * qsplitrainparamb * qsplitrainparamc * (one - q(:)) / &
1073  ((qsplitrainparamc * qnv(:)) + one) ** 2
1074  qi(:) = (one - lf(:)) * qsplitrainparama * qsplitrainparamb * qsplitrainparamc * (one - q(:)) / &
1075  ((qsplitrainparamc * qnv(:)) + one) ** 2
1076 
1077  else
1078 
1079  ql(:) = lf(:) * (one - q(:))
1080  qi(:) = (one - lf(:)) * (one - q(:))
1081 
1082  end if
1083 
1084  end if
1085 
1086 end subroutine ops_satrad_qsplit
1087 
1088 !-------------------------------------------------------------------------------
1089 !> Do cholesky decomposition
1090 !!
1091 !! \details Heritage: Ops_Cholesky.inc
1092 !!
1093 !! Solves the Linear equation UQ=V for Q where U is a symmetric positive definite
1094 !! matrix and U and Q are vectors of length N. The method follows that in Golub
1095 !! and Van Loan although this is pretty standard.
1096 !!
1097 !! if U is not positive definite this will be detected by the program and flagged
1098 !! as an error. U is assumed to be symmetric as only the upper triangle is in
1099 !! fact used.
1100 !!
1101 !! \author Met Office
1102 !!
1103 !! \date 09/06/2020: Created
1104 !!
1105 subroutine ops_cholesky (U, &
1106  V, &
1107  N, &
1108  Q, &
1109  ErrorCode)
1110 
1111 implicit none
1112 
1113 ! subroutine arguments:
1114 integer, intent(in) :: n
1115 real(kind_real), intent(in) :: u(n,n)
1116 real(kind_real), intent(in) :: v(n)
1117 real(kind_real), intent(out) :: q(n)
1118 integer, intent(out) :: errorcode
1119 
1120 ! Local declarations:
1121 character(len=*), parameter :: routinename = "Ops_Cholesky"
1122 real(kind_real), parameter :: tolerance = tiny (zero) * 100.0_kind_real
1123 character(len=50) :: errormessage
1124 integer :: j
1125 integer :: k
1126 real(kind_real) :: g(n,n) ! The Cholesky Triangle Matrix
1127 real(kind_real) :: x(n) ! Temporary array used in calculating G
1128 
1129 errorcode = 0
1130 
1131 ! Determine the Cholesky triangle matrix.
1132 
1133 do j = 1, n
1134  x(j:n) = u(j:n,j)
1135  if (j /= 1) then
1136  do k = 1, j - 1
1137  x(j:n) = x(j:n) - g(j,k) * g(j:n,k)
1138  end do
1139  end if
1140  if (x(j) <= tolerance) then
1141  errorcode = 1
1142  errormessage = ' :U matrix is not positive definite'
1143  write(*,*) trim(routinename), trim(errormessage)
1144  goto 9999
1145  end if
1146  g(j:n,j) = x(j:n) / sqrt(x(j))
1147 end do
1148 
1149 ! Solve Gx=v for x by forward substitution
1150 
1151 x = v
1152 x(1) = x(1) / g(1,1)
1153 do j = 2, n
1154  x(j) = (x(j) - dot_product(g(j,1:j - 1), x(1:j - 1))) / g(j,j)
1155 end do
1156 
1157 ! Solve G^T.q=x for q by backward substitution
1158 
1159 q = x
1160 q(n) = q(n) / g(n,n)
1161 do j = n - 1, 1, -1
1162  q(j) = (q(j) - dot_product(g(j + 1:n,j), q(j + 1:n))) / g(j,j)
1163 end do
1164 
1165 9999 continue
1166 
1167 end subroutine ops_cholesky
1168 
1169 !-------------------------------------------------------------------------------
1170 !> Find a free file unit.
1171 !!
1172 !! \author Met Office
1173 !!
1174 !! \date 09/06/2020: Created
1175 !!
1176 function ufo_utils_iogetfreeunit() result(unit)
1177 
1178 implicit none
1179 
1180 integer :: unit
1181 
1182 integer, parameter :: unit_min=10
1183 integer, parameter :: unit_max=1000
1184 logical :: opened
1185 integer :: lun
1186 integer :: newunit
1187 
1188 newunit=-1
1189 do lun=unit_min,unit_max
1190  inquire(unit=lun,opened=opened)
1191  if (.not. opened) then
1192  newunit=lun
1193  exit
1194  end if
1195 end do
1196 unit=newunit
1197 
1198 end function ufo_utils_iogetfreeunit
1199 
1200 !-------------------------------------------------------------------------------
1201 ! Invert a matrix.
1202 !
1203 ! Variables with intent in:
1204 !
1205 ! N: Size of the matrix being inverted
1206 ! M: If MATRIX is not present this is the same as N, else this is the other
1207 ! dimension of MATRIX.
1208 !
1209 ! Variables with intent inout:
1210 !
1211 ! A: Real matrix (assumed square and symmetrical) overwritten by its
1212 ! inverse if MATRIX is not present.
1213 !
1214 ! Variables with intent out:
1215 !
1216 ! ExitCode: 0: ok, 1: A is not positive definite.
1217 !
1218 ! Optional variables with intent inout:
1219 !
1220 ! Matrix: If present this input matrix is replaced by (Matrix).A^-1 on exit
1221 ! (leaving A unchanged).
1222 !
1223 ! Uses Cholesky decomposition - a method particularly suitable for real
1224 ! symmetric matrices. Cholesky decomposition solves the Linear equation UQ=V
1225 ! for Q where U is a symmetric positive definite matrix and U and Q are vectors
1226 ! of length N. The method follows that in Golub and Van Loan although this is
1227 ! pretty standard. If U is not positive definite this will be detected by the
1228 ! program and flagged as an error. U is assumed to be symmetric as only the
1229 ! upper triangle is in fact used. If the the optional parameter Matrix is
1230 ! present, it is replaced by (Matrix).A^-1 on exit and A is left unchanged.
1231 !-------------------------------------------------------------------------------
1232 
1233 subroutine invertmatrix (n, &
1234  m, &
1235  a, &
1236  status, &
1237  matrix)
1238 
1239 implicit none
1240 
1241 ! subroutine arguments:
1242 integer, intent(in) :: n !< order of a
1243 integer, intent(in) :: m !< order of optional matrix, if required
1244 real(kind=kind_real), intent(inout) :: a(n,n) !< square mx, overwritten by its inverse
1245 integer, intent(out) :: status !< 0 if all ok 1 if matrix is not positive definite
1246 real(kind=kind_real), optional, intent(inout) :: matrix(n,m) !< replaced by (matrix).a^-1 on exit
1247 
1248 ! local declarations:
1249 character(len=*), parameter :: routinename = "InvertMatrix"
1250 character(len=80) :: errormessage
1251 real(kind=kind_real), parameter :: tolerance = tiny(zero) * 100.0_kind_real
1252 integer :: i
1253 integer :: j
1254 integer :: k
1255 integer :: mm
1256 real(kind=kind_real) :: g(n,n) ! the cholesky triangle matrix
1257 real(kind=kind_real) :: q(n)
1258 real(kind=kind_real) :: tmp(n,m)
1259 real(kind=kind_real) :: v(n)
1260 real(kind=kind_real) :: x(n) ! temporary array
1261 
1262 status = 0
1263 
1264 ! determine the cholesky triangle matrix.
1265 
1266 do j = 1, n
1267  x(j:n) = a(j:n,j)
1268  if (j /= 1) then
1269  do k = 1, j - 1
1270  x(j:n) = x(j:n) - g(j,k) * g(j:n,k)
1271  end do
1272  end if
1273  if (x(j) <= tolerance) then
1274  errormessage = 'Matrix is not positive definite'
1275  call fckit_log % warning(errormessage)
1276  status = 1
1277  goto 9999
1278  end if
1279  g(j:n,j) = x(j:n) / sqrt(x(j))
1280 end do
1281 
1282 ! now solve the equation g.g^t.q=v for the set of
1283 ! vectors, v, with one element = 1 and the rest zero.
1284 ! the solutions q are brought together at the end to form
1285 ! the inverse of g.g^t (i.e., the inverse of a).
1286 
1287 if (present (matrix)) then
1288  mm = m
1289 else
1290  ! make sure that the dimensions of tmp were correctly specified
1291  if (m /= n) then
1292  errormessage = '2nd and 3rd arguments of routine must be'
1293  call fckit_log % warning(errormessage)
1294  errormessage = 'identical if the matrix option is not present'
1295  call fckit_log % warning(errormessage)
1296  status = 2
1297  goto 9999
1298  end if
1299  mm = n
1300 end if
1301 
1302 main_loop : do i = 1, mm
1303  if (.not. (present (matrix))) then
1304  v(:) = zero
1305  v(i) = one
1306  else
1307  v(:) = matrix(i,:)
1308  end if
1309 
1310  ! solve gx=v for x by forward substitution
1311 
1312  x = v
1313  x(1) = x(1) / g(1,1)
1314  do j = 2, n
1315  x(j) = (x(j) - dot_product(g(j,1:j - 1), x(1:j - 1))) / g(j,j)
1316  end do
1317 
1318  ! solve g^t.q=x for q by backward substitution
1319 
1320  q = x
1321  q(n) = q(n) / g(n,n)
1322  do j = n - 1, 1, -1
1323  q(j) = (q(j) - dot_product(g(j + 1:n,j), q(j + 1:n))) / g(j,j)
1324  end do
1325 
1326  tmp(:,i) = q(:)
1327 
1328 end do main_loop
1329 
1330 if (.not. (present (matrix))) then
1331  a(:,:) = tmp(:,:)
1332 else
1333  matrix(:,:) = tmp(:,:)
1334 end if
1335 
1336 9999 continue
1337 
1338 end subroutine invertmatrix
1339 
1340 !-------------------------------------------------------------------------------
1341 ! This function will comapre two strings while avoiding trim in order to speed
1342 ! up the string comparison. This will allow for trailing blanks on either string
1343 ! to count for a match.
1344 !
1345 ! Since the strings may be different lengths coming in (due to trailing blanks)
1346 ! a helper function (cmp_ordered_strings) is called with the shorter string first
1347 ! which helps simplify the logic of the comparison.
1348 
1349 function cmp_strings(str1, str2)
1350  implicit none
1351 
1352  logical :: cmp_strings
1353  character(len=*) :: str1
1354  character(len=*) :: str2
1355 
1356  if (len(str1) <= len(str2)) then
1357  cmp_strings = cmp_ordered_strings(str1, str2)
1358  else
1359  cmp_strings = cmp_ordered_strings(str2, str1)
1360  endif
1361 
1362  return
1363 end function cmp_strings
1364 
1365 !-------------------------------------------------------------------------------
1366 ! This is a helper function for the public cmp_strings function. It is intended
1367 ! for this function to be called with the arguments of cmp_strings, with the shorter
1368 ! of the two strings in the first argument.
1369 !
1370 ! The algorithm is to compare each character one by one between each string. Then if
1371 ! one string is longer, continue looking at that string an make sure the remainder
1372 ! consists of blanks before calling it a match between the strings. Knowing that the
1373 ! first argument is the shorter of the two stings helps simplify the code doing the check.
1374 !
1375 ! The trim call is avoided since it allocates, copies strings and deallocates which
1376 ! collectively are operations that are too expensive.
1377 
1378 function cmp_ordered_strings(shorter_str, longer_str)
1379  implicit none
1380 
1381  logical :: cmp_ordered_strings
1382  character(len=*) :: shorter_str
1383  character(len=*) :: longer_str
1384 
1385  integer :: i
1386  integer :: j
1387 
1388  ! Check each character from the beginning of the strings to the length of the
1389  ! shorter string. Exit the loop right away if a mis-match occurs to save time.
1390  cmp_ordered_strings = .true.
1391  do i = 1, len(shorter_str)
1392  if (shorter_str(i:i) /= longer_str(i:i)) then
1393  cmp_ordered_strings = .false.
1394  exit
1395  endif
1396  enddo
1397 
1398  ! If cmp_ordered_strings is false we can return immediately. If true, then we need to check
1399  ! that the remainder of the longer string constists of all blank spaces before
1400  ! declaring the strings equal.
1401  if (cmp_ordered_strings) then
1402  do j = i, len(longer_str)
1403  if (longer_str(j:j) /= " ") then
1404  cmp_ordered_strings = .false.
1405  exit
1406  endif
1407  enddo
1408  endif
1409 
1410  return
1411 end function cmp_ordered_strings
1412 
1413 !-------------------------------------------------------------------------------
1414 
1415 end module ufo_utils_mod
real(kind_real), parameter, public min_q
real(kind_real), parameter, public one
real(kind_real), parameter, public t0c
real(kind_real), parameter, public zero
real(kind_real), parameter, public zerodegc
real(kind_real), parameter, public two
real(kind_real), parameter, public half
real(kind_real), parameter, public epsilon
Fortran module with various useful routines.
logical function cmp_ordered_strings(shorter_str, longer_str)
subroutine, public ops_satrad_qsplit(output_type, p, t, qtotal, q, ql, qi, UseQtSplitRain)
Split the humidity into water vapour, liquid water and ice.
subroutine, public ops_qsat(QS, T, P, npnts)
Calculate the Saturation Specific Humidity Scheme (Qsat): Vapour to Liquid/Ice.
subroutine, public ops_qsatwat(QS, T, P, npnts)
Saturation Specific Humidity Scheme: Vapour to Liquid.
integer function, public ufo_utils_iogetfreeunit()
Find a free file unit.
logical function, public cmp_strings(str1, str2)
subroutine, public ops_cholesky(U, V, N, Q, ErrorCode)
Do cholesky decomposition.
subroutine, public invertmatrix(n, m, a, status, matrix)