1module humidityLimits_mod
2 ! MODULE humidityLimits_mod (prefix='qlim' category='4. Data Object transformations')
3 !
4 !:Purpose: Various manipulations of humidity-related quantities.
5 !
6 use midasMpi_mod
7 use utilities_mod
8 use mathPhysConstants_mod
9 use varNameList_mod
10 use physicsFunctions_mod
11 use verticalCoord_mod
12 use gridStateVector_mod
13 use ensembleStateVector_mod
14 use calcHeightAndPressure_mod
15 implicit none
16 save
17 private
18
19 ! public procedures
20 public :: qlim_saturationLimit, qlim_rttovLimit, qlim_setMin
21 public :: qlim_getMinValueCloud, qlim_getMaxValueCloud
22
23 real(8), parameter :: mixratio_to_ppmv = 1.60771704d+6
24 real(8) :: qlim_minValueLWCR, qlim_minValueIWCR, qlim_minValueRF, qlim_minValueSF, qlim_minValueCLDR
25 real(8) :: qlim_maxValueLWCR, qlim_maxValueIWCR, qlim_maxValueRF, qlim_maxValueSF, qlim_maxValueCLDR
26
27 ! interface for qlim_saturationLimit
28 interface qlim_saturationLimit
29 module procedure qlim_saturationLimit_gsv
30 module procedure qlim_saturationLimit_ens
31 end interface qlim_saturationLimit
32
33 ! interface for qlim_rttovLimit
34 interface qlim_rttovLimit
35 module procedure qlim_rttovLimit_gsv
36 module procedure qlim_rttovLimit_ens
37 end interface qlim_rttovLimit
38
39 ! interface for qlim_setMin
40 interface qlim_setMin
41 module procedure qlim_setMin_ens
42 end interface qlim_setMin
43
44contains
45
46 !--------------------------------------------------------------------------
47 ! readNameList
48 !--------------------------------------------------------------------------
49 subroutine readNameList
50 !
51 !:Purpose: Reading NAMQLIM namelist by any subroutines in humidityLimits_mod module.
52 !
53 implicit none
54
55 ! Locals:
56 integer :: nulnam, ierr
57 integer, external :: fnom, fclos
58 logical, save :: nmlAlreadyRead = .false.
59
60 ! Namelist variables:
61 real(8) :: minValueLWCR ! minimum LWCR value
62 real(8) :: maxValueLWCR ! maximum LWCR value
63 real(8) :: minValueIWCR ! minimum IWCR value
64 real(8) :: maxValueIWCR ! maximum IWCR value
65 real(8) :: minValueRF ! minimum RF value
66 real(8) :: maxValueRF ! maximum RF value
67 real(8) :: minValueSF ! minimum SF value
68 real(8) :: maxValueSF ! maximum SF value
69 real(8) :: minValueCLDR ! minimum CLDR value
70 real(8) :: maxValueCLDR ! maximum CLDR value
71
72 NAMELIST /NAMQLIM/ minValueLWCR, maxValueLWCR, minValueIWCR, maxValueIWCR, &
73 minValueRF, maxValueRF, minValueSF, maxValueSF, minValueCLDR, maxValueCLDR
74
75 if ( nmlAlreadyRead ) return
76
77 nmlAlreadyRead = .true.
78
79 !- Setting default values
80 minValueLWCR = 1.0d-9
81 maxValueLWCR = 1.0d0
82
83 minValueIWCR = 1.0d-9
84 maxValueIWCR = 1.0d0
85
86 minValueRF = 0.0d0
87 maxValueRF = 1.0d0
88
89 minValueSF = 0.0d0
90 maxValueSF = 1.0d0
91
92 minValueCLDR = 0.0d0
93 maxValueCLDR = 1.0d0
94
95 if ( .not. utl_isNamelistPresent('NAMQLIM','./flnml') ) then
96 if ( mmpi_myid == 0 ) then
97 write(*,*) 'NAMQLIM is missing in the namelist. The default values will be taken.'
98 end if
99
100 else
101 ! Reading the namelist
102 nulnam = 0
103 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
104 read(nulnam, nml=namqlim, iostat=ierr)
105 if ( ierr /= 0) call utl_abort('humidityLimits_mod: Error reading namelist')
106 ierr = fclos(nulnam)
107
108 end if
109 if ( mmpi_myid == 0 ) write(*,nml=namqlim)
110
111 ! Transfer namelist variables to module variables.
112 qlim_minValueLWCR = minValueLWCR
113 qlim_maxValueLWCR = maxValueLWCR
114
115 qlim_minValueIWCR = minValueIWCR
116 qlim_maxValueIWCR = maxValueIWCR
117
118 qlim_minValueRF = minValueRF
119 qlim_maxValueRF = maxValueRF
120
121 qlim_minValueSF = minValueSF
122 qlim_maxValueSF = maxValueSF
123
124 qlim_minValueCLDR = minValueCLDR
125 qlim_maxValueCLDR = maxValueCLDR
126
127 end subroutine readNameList
128
129 !--------------------------------------------------------------------------
130 ! qlim_saturationLimit_gsv
131 !--------------------------------------------------------------------------
132 subroutine qlim_saturationLimit_gsv(statevector)
133 !
134 !:Purpose: To impose saturation limit on humidity variable of a statevector
135 !
136 implicit none
137
138 ! Arguments:
139 type(struct_gsv), intent(inout) :: statevector
140
141 ! Locals:
142 type(struct_vco), pointer :: vco_ptr
143 real(4), pointer :: hu_ptr_r4(:,:,:,:), tt_ptr_r4(:,:,:,:), psfc_ptr_r4(:,:,:,:), psfcLS_ptr_r4(:,:,:,:)
144 real(8), pointer :: hu_ptr_r8(:,:,:,:), tt_ptr_r8(:,:,:,:), psfc_ptr_r8(:,:,:,:), psfcLS_ptr_r8(:,:,:,:)
145 real(8), pointer :: pressure(:,:,:)
146 real(8) :: hu, husat, hu_modified, tt
147 integer :: lon1, lon2, lat1, lat2, lev1, lev2
148 real(8), allocatable :: psfc(:,:), psfcLS(:,:)
149 integer :: lonIndex, latIndex, levIndex, stepIndex
150
151 if (mmpi_myid == 0) write(*,*) 'qlim_saturationLimit_gsv: STARTING'
152
153 if( .not. gsv_varExist(statevector,'HU') ) then
154 if( mmpi_myid == 0 ) write(*,*) 'qlim_saturationLimit_gsv: statevector does not ' // &
155 'contain humidity ... doing nothing'
156 return
157 end if
158
159 vco_ptr => gsv_getVco(statevector)
160 if (stateVector%dataKind == 8) then
161 call gsv_getField(statevector,hu_ptr_r8,'HU')
162 call gsv_getField(statevector,tt_ptr_r8,'TT')
163 else
164 call gsv_getField(statevector,hu_ptr_r4,'HU')
165 call gsv_getField(statevector,tt_ptr_r4,'TT')
166 end if
167
168 lon1 = statevector%myLonBeg
169 lon2 = statevector%myLonEnd
170 lat1 = statevector%myLatBeg
171 lat2 = statevector%myLatEnd
172 lev1 = 1
173 lev2 = gsv_getNumLev(statevector,'TH')
174
175 allocate(psfc(lon2-lon1+1,lat2-lat1+1))
176 if (vco_ptr%vcode == 5100) allocate(psfcLS(lon2-lon1+1,lat2-lat1+1))
177 do stepIndex = 1, statevector%numStep
178 if (stateVector%dataKind == 8) then
179 call gsv_getField(statevector,psfc_ptr_r8,'P0')
180 psfc(:,:) = psfc_ptr_r8(:,:,1,stepIndex)
181 if (vco_ptr%vcode == 5100) then
182 call gsv_getField(statevector,psfcLS_ptr_r8,'P0LS')
183 psfcLS(:,:) = psfcLS_ptr_r8(:,:,1,stepIndex)
184 end if
185 else
186 call gsv_getField(statevector,psfc_ptr_r4,'P0')
187 psfc(:,:) = psfc_ptr_r4(:,:,1,stepIndex)
188 if (vco_ptr%vcode == 5100) then
189 call gsv_getField(statevector,psfcLS_ptr_r4,'P0LS')
190 psfcLS(:,:) = psfcLS_ptr_r4(:,:,1,stepIndex)
191 end if
192 end if
193 if (vco_ptr%vcode == 5100) then
194 call czp_fetch3DLevels(vco_ptr, psfc, sfcFldLS_opt=psfcLS, fldT_opt=pressure)
195 else
196 call czp_fetch3DLevels(vco_ptr, psfc, fldT_opt=pressure)
197 end if
198 if (stateVector%dataKind == 8) then
199 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex, hu, tt, husat, hu_modified)
200 do levIndex = lev1, lev2
201 do latIndex = lat1, lat2
202 do lonIndex = lon1, lon2
203 hu = hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex)
204 tt = tt_ptr_r8(lonIndex,latIndex,levIndex,stepIndex)
205
206 ! get the saturated vapor pressure from HU
207 husat = phf_foqst8(tt, pressure(lonIndex-lon1+1,latIndex-lat1+1,levIndex) )
208
209 ! limit the humidity to the saturated humidity
210 hu_modified = min(husat, hu)
211 hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = hu_modified
212
213 end do ! lonIndex
214 end do ! latIndex
215 end do ! levIndex
216 !$OMP END PARALLEL DO
217 else
218 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex, hu, tt, husat, hu_modified)
219 do levIndex = lev1, lev2
220 do latIndex = lat1, lat2
221 do lonIndex = lon1, lon2
222 hu = hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex)
223 tt = tt_ptr_r4(lonIndex,latIndex,levIndex,stepIndex)
224
225 ! get the saturated vapor pressure from HU
226 husat = phf_foqst8(tt, pressure(lonIndex-lon1+1,latIndex-lat1+1,levIndex) )
227
228 ! limit the humidity to the saturated humidity
229 hu_modified = min(husat, hu)
230 hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = hu_modified
231
232 end do ! lonIndex
233 end do ! latIndex
234 end do ! levIndex
235 !$OMP END PARALLEL DO
236 end if
237
238 deallocate(pressure)
239
240 end do ! stepIndex
241
242 deallocate(psfc)
243 if (allocated(psfcLS)) deallocate(psfcLS)
244
245 end subroutine qlim_saturationLimit_gsv
246
247 !--------------------------------------------------------------------------
248 ! qlim_saturationLimit_ens
249 !--------------------------------------------------------------------------
250 subroutine qlim_saturationLimit_ens(ensemble)
251 !
252 !:Purpose: To impose saturation limit on humidity variable of an ensemble
253 !
254 implicit none
255
256 ! Arguments:
257 type(struct_ens), intent(inout) :: ensemble
258
259 ! Locals:
260 type(struct_vco), pointer :: vco_ptr
261 real(4), pointer :: hu_ptr_r4(:,:,:,:), tt_ptr_r4(:,:,:,:), psfc_ptr_r4(:,:,:,:), psfcLS_ptr_r4(:,:,:,:)
262 real(8), pointer :: pressure(:,:,:)
263 real(8) :: hu, husat, hu_modified, tt
264 integer :: lon1, lon2, lat1, lat2, numLev
265 real(8), allocatable :: psfc(:,:),psfcLS(:,:)
266 integer :: lonIndex, latIndex, levIndex, stepIndex, memberIndex, varLevIndex
267
268 if (mmpi_myid == 0) write(*,*) 'qlim_saturationLimit_ens: STARTING'
269
270 if (ens_getDataKind(ensemble) == 8) then
271 call utl_abort('qlim_saturationLimit_ens: Not compatible with dataKind = 8')
272 end if
273
274 if( .not. ens_varExist(ensemble,'HU') ) then
275 if( mmpi_myid == 0 ) write(*,*) 'qlim_saturationLimit_ens: ensemble does not ' // &
276 'contain humidity ... doing nothing'
277 return
278 end if
279
280 vco_ptr => ens_getVco(ensemble)
281 numLev = ens_getNumLev(ensemble,'TH')
282 call ens_getLatLonBounds(ensemble, lon1, lon2, lat1, lat2)
283 allocate(psfc(ens_getNumMembers(ensemble),ens_getNumStep(ensemble)))
284 if (vco_ptr%vcode == 5100) allocate(psfcLS(ens_getNumMembers(ensemble),ens_getNumStep(ensemble)))
285
286 do latIndex = lat1, lat2
287 do lonIndex = lon1, lon2
288
289 ! compute pressure for all members and steps
290 varLevIndex = ens_getKFromLevVarName(ensemble, 1, 'P0')
291 psfc_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
292 psfc(:,:) = psfc_ptr_r4(:,:,lonIndex,latIndex)
293 if (vco_ptr%vcode == 5100) then
294 varLevIndex = ens_getKFromLevVarName(ensemble, 1, 'P0LS')
295 psfcLS_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
296 psfcLS(:,:) = psfcLS_ptr_r4(:,:,lonIndex,latIndex)
297 call czp_fetch3DLevels(vco_ptr, psfc, sfcFldLS_opt=psfcLS, fldT_opt=pressure)
298 else
299 call czp_fetch3DLevels(vco_ptr, psfc, fldT_opt=pressure)
300 end if
301
302 do levIndex = 1, numLev
303 varLevIndex = ens_getKFromLevVarName(ensemble, levIndex, 'HU')
304 hu_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
305 varLevIndex = ens_getKFromLevVarName(ensemble, levIndex, 'TT')
306 tt_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
307
308 !$OMP PARALLEL DO PRIVATE (stepIndex, memberIndex, hu, tt, husat, hu_modified)
309 do stepIndex = 1, ens_getNumStep(ensemble)
310 do memberIndex = 1, ens_getNumMembers(ensemble)
311 hu = hu_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex)
312 tt = tt_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex)
313
314 ! get the saturated vapor pressure from HU
315 husat = phf_foqst8(tt, pressure(memberIndex,stepIndex,levIndex) )
316
317 ! limit the humidity to the saturated humidity
318 hu_modified = min(husat, hu)
319 hu_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = hu_modified
320
321 end do ! memberIndex
322 end do ! stepIndex
323 !$OMP END PARALLEL DO
324
325 end do ! levIndex
326 deallocate(pressure)
327
328 end do ! lonIndex
329 end do ! latIndex
330
331 deallocate(psfc)
332 if (allocated(psfcLS)) deallocate(psfcLS)
333
334 end subroutine qlim_saturationLimit_ens
335
336 !--------------------------------------------------------------------------
337 ! qlim_rttovLimit_gsv
338 !--------------------------------------------------------------------------
339 subroutine qlim_rttovLimit_gsv(statevector, varName_opt, applyLimitToCloud_opt)
340 !
341 !:Purpose: To impose RTTOV limits on humidity/cloud
342 !
343 implicit none
344
345 ! Arguments:
346 type(struct_gsv), intent(inout) :: statevector
347 character(len=*), optional, intent(in) :: varName_opt
348 logical, optional, intent(in) :: applyLimitToCloud_opt
349
350 ! Locals:
351 type(struct_vco), pointer :: vco_ptr
352 real(8), allocatable :: press_rttov(:), qmin_rttov(:), qmax_rttov(:)
353 real(8), allocatable :: psfc(:,:),psfcLS(:,:)
354 real(8), allocatable :: qmin3D_rttov(:,:,:), qmax3D_rttov(:,:,:)
355 real(8), pointer :: hu_ptr_r8(:,:,:,:), psfc_ptr_r8(:,:,:,:), psfcLS_ptr_r8(:,:,:,:)
356 real(8), pointer :: cld_ptr_r8(:,:,:,:)
357 real(4), pointer :: hu_ptr_r4(:,:,:,:), psfc_ptr_r4(:,:,:,:), psfcLS_ptr_r4(:,:,:,:)
358 real(4), pointer :: cld_ptr_r4(:,:,:,:)
359 real(8), pointer :: pressure(:,:,:)
360 real(8) :: hu, hu_modified
361 real(8) :: cld, cld_modified
362 real(8) :: minValueCld, maxValueCld
363 integer :: lon1, lon2, lat1, lat2, lev1, lev2, varNameIndex
364 integer :: lonIndex, latIndex, levIndex, stepIndex
365 integer :: ni, nj, numLev, numLev_rttov
366 integer :: fnom, fclos, ierr, nulfile
367 character(len=256) :: fileName
368 character(len=4) :: varName
369 logical, save :: firstTime=.true.
370 logical :: applyLimitToAllVarname
371 logical :: applyLimitToCloud
372
373 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_gsv: STARTING'
374
375 if (present(varName_opt)) then
376 applyLimitToAllVarname = .false.
377 varName = varName_opt
378 else
379 applyLimitToAllVarname = .true.
380 varName = 'XXXX'
381 end if
382
383 ! for cloud, limits are applied to ALL cloud variables.
384 if (present(applyLimitToCloud_opt)) then
385 applyLimitToCloud = applyLimitToCloud_opt
386 if (.not. applyLimitToCloud) call utl_abort('qlim_rttovLimit_gsv: remove applyLimitToCloud_opt argument')
387 else
388 if (vnl_isCloudVar(varName)) then
389 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_gsv: limits are applied to ALL cloud variables'
390 applyLimitToCloud = .true.
391 else
392 applyLimitToCloud = .false.
393 end if
394 end if
395
396 if (applyLimitToCloud) applyLimitToAllVarname = .false.
397
398 if ((applyLimitToAllVarname .or. trim(varName) == 'HU') .and. &
399 gsv_varExist(statevector,'HU')) then
400
401 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_gsv: applying limits to HU.'
402
403 ! Read in RTTOV humidity limits
404 fileName = "rttov_h2o_limits.dat"
405 nulfile = 0
406 ierr = fnom(nulfile, fileName, "FMT+OLD+R/O", 0)
407 if( ierr /= 0 ) then
408 if ( mmpi_myid == 0 ) write(*,*) 'fileName = ', fileName
409 call utl_abort('qlim_rttovLimit_gsv: error opening the humidity limits file')
410 end if
411
412 read(nulfile,*) numLev_rttov
413 if ( mmpi_myid == 0 .and. firstTime ) write(*,*) 'qlim_rttovLimit_gsv: rttov number of levels = ', numLev_rttov
414 allocate(press_rttov(numLev_rttov))
415 allocate(qmin_rttov(numLev_rttov))
416 allocate(qmax_rttov(numLev_rttov))
417 do levIndex = 1, numLev_rttov
418 read(nulfile,*) press_rttov(levIndex), qmax_rttov(levIndex), qmin_rttov(levIndex)
419 end do
420 ierr = fclos(nulfile)
421 press_rttov(:) = press_rttov(:) * mpc_pa_per_mbar_r8
422 qmin_rttov(:) = qmin_rttov(:) / mixratio_to_ppmv
423 qmax_rttov(:) = qmax_rttov(:) / mixratio_to_ppmv
424
425 if (firstTime) then
426 write(*,*) ' '
427 do levIndex = 1, numLev_rttov
428 if ( mmpi_myid == 0 ) write(*,fmt='(" qlim_rttovLimit_gsv: LEVEL = ",I4,", PRES = ",F9.0,", HUMIN = ",E10.2,", HUMAX = ",E10.2)') &
429 levIndex, press_rttov(levIndex), qmin_rttov(levIndex), qmax_rttov(levIndex)
430 end do
431 firstTime = .false.
432 end if
433
434 vco_ptr => gsv_getVco(statevector)
435 if (statevector%dataKind == 8) then
436 call gsv_getField(statevector,hu_ptr_r8,'HU')
437 else
438 call gsv_getField(statevector,hu_ptr_r4,'HU')
439 end if
440
441 lon1 = statevector%myLonBeg
442 lon2 = statevector%myLonEnd
443 lat1 = statevector%myLatBeg
444 lat2 = statevector%myLatEnd
445 lev1 = 1
446 lev2 = gsv_getNumLev(statevector,'TH')
447
448 ni = lon2 - lon1 + 1
449 nj = lat2 - lat1 + 1
450 numLev = lev2 - lev1 + 1
451 allocate( qmin3D_rttov(ni,nj,numLev) )
452 allocate( qmax3D_rttov(ni,nj,numLev) )
453 allocate( psfc(lon2-lon1+1,lat2-lat1+1) )
454 if (vco_ptr%vcode == 5100) allocate( psfcLS(lon2-lon1+1,lat2-lat1+1) )
455
456 do stepIndex = 1, statevector%numStep
457 if (statevector%dataKind == 8) then
458 call gsv_getField(statevector,psfc_ptr_r8,'P0')
459 psfc(:,:) = psfc_ptr_r8(:,:,1,stepIndex)
460 if (vco_ptr%vcode == 5100) then
461 call gsv_getField(statevector,psfcLS_ptr_r8,'P0')
462 psfcLS(:,:) = psfcLS_ptr_r8(:,:,1,stepIndex)
463 end if
464 else
465 call gsv_getField(statevector,psfc_ptr_r4,'P0')
466 psfc(:,:) = real(psfc_ptr_r4(:,:,1,stepIndex),8)
467 if (vco_ptr%vcode == 5100) then
468 call gsv_getField(statevector,psfcLS_ptr_r4,'P0')
469 psfcLS(:,:) = psfcLS_ptr_r4(:,:,1,stepIndex)
470 end if
471 end if
472 if (vco_ptr%vcode == 5100) then
473 call czp_fetch3DLevels(vco_ptr, psfc, sfcFldLS_opt=psfcLS, fldT_opt=pressure)
474 else
475 call czp_fetch3DLevels(vco_ptr, psfc, fldT_opt=pressure)
476 end if
477
478 ! Interpolate RTTOV limits onto model levels
479 call qlim_lintv_minmax(press_rttov, qmin_rttov, qmax_rttov, numLev_rttov, &
480 ni, nj, numLev, pressure, qmin3D_rttov, qmax3D_rttov)
481
482 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex, hu, hu_modified)
483 do levIndex = lev1, lev2
484 do latIndex = lat1, lat2
485 do lonIndex = lon1, lon2
486 if (statevector%dataKind == 8) then
487 hu = hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex)
488 else
489 hu = real(hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex),8)
490 end if
491
492 ! limit the humidity according to the rttov limits
493 hu_modified = max(hu, qmin3D_rttov(lonIndex - lon1 + 1, latIndex - lat1 + 1, levIndex) )
494 hu_modified = min(hu_modified, qmax3D_rttov(lonIndex - lon1 + 1, latIndex - lat1 + 1, levIndex) )
495 if (statevector%dataKind == 8) then
496 hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = hu_modified
497 else
498 hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = real(hu_modified,4)
499 end if
500
501 end do ! lonIndex
502 end do ! latIndex
503 end do ! levIndex
504 !$OMP END PARALLEL DO
505
506 deallocate(pressure)
507 end do ! stepIndex
508
509 deallocate( psfc )
510 if (allocated(psfcLS)) deallocate( psfcLS )
511 deallocate( qmin3D_rttov )
512 deallocate( qmax3D_rttov )
513
514 deallocate(qmax_rttov)
515 deallocate(qmin_rttov)
516 deallocate(press_rttov)
517
518 end if
519
520 ! apply limits to ALL available cloud variables
521 if ((applyLimitToAllVarname .or. applyLimitToCloud) .and. cloudExistInStateVector(statevector)) then
522
523 do varNameIndex = 1, vnl_numvarmaxCloud
524 if (.not. gsv_varExist(statevector, vnl_varNameListCloud(varNameIndex))) cycle
525
526 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_gsv: applying limits to ', &
527 vnl_varNameListCloud(varNameIndex)
528
529 if (statevector%dataKind == 8) then
530 call gsv_getField(statevector, cld_ptr_r8, vnl_varNameListCloud(varNameIndex))
531 else
532 call gsv_getField(statevector, cld_ptr_r4, vnl_varNameListCloud(varNameIndex))
533 end if
534
535 minValueCld = qlim_getMinValueCloud(vnl_varNameListCloud(varNameIndex))
536 maxValueCld = qlim_getMaxValueCloud(vnl_varNameListCloud(varNameIndex))
537
538 lon1 = statevector%myLonBeg
539 lon2 = statevector%myLonEnd
540 lat1 = statevector%myLatBeg
541 lat2 = statevector%myLatEnd
542 lev1 = 1
543 lev2 = gsv_getNumLev(statevector,'TH')
544
545 do stepIndex = 1, statevector%numStep
546
547 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex, cld, cld_modified)
548 do levIndex = lev1, lev2
549 do latIndex = lat1, lat2
550 do lonIndex = lon1, lon2
551 if (statevector%dataKind == 8) then
552 cld = cld_ptr_r8(lonIndex,latIndex,levIndex,stepIndex)
553 else
554 cld = real(cld_ptr_r4(lonIndex,latIndex,levIndex,stepIndex),8)
555 end if
556
557 cld_modified = max(cld,minValueCld)
558 cld_modified = min(cld_modified,maxValueCld)
559 if (statevector%dataKind == 8) then
560 cld_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = cld_modified
561 else
562 cld_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = real(cld_modified,4)
563 end if
564
565 end do ! lonIndex
566 end do ! latIndex
567 end do ! levIndex
568 !$OMP END PARALLEL DO
569
570 end do ! stepIndex
571 end do ! varNameIndex
572
573 end if
574
575 end subroutine qlim_rttovLimit_gsv
576
577 !--------------------------------------------------------------------------
578 ! qlim_rttovLimit_ens
579 !--------------------------------------------------------------------------
580 subroutine qlim_rttovLimit_ens(ensemble, varName_opt, applyLimitToCloud_opt)
581 !
582 !:Purpose: To impose RTTOV limits on humidity/cloud
583 !
584 implicit none
585
586 ! Arguments:
587 type(struct_ens), intent(inout) :: ensemble
588 character(len=*), optional, intent(in) :: varName_opt
589 logical, optional, intent(in) :: applyLimitToCloud_opt
590
591 ! Locals:
592 type(struct_vco), pointer :: vco_ptr
593 real(8), allocatable :: press_rttov(:), qmin_rttov(:), qmax_rttov(:)
594 real(8), allocatable :: psfc(:,:),psfcLS(:,:)
595 real(8), allocatable :: qmin3D_rttov(:,:,:), qmax3D_rttov(:,:,:)
596 real(4), pointer :: hu_ptr_r4(:,:,:,:), psfc_ptr_r4(:,:,:,:), psfcLS_ptr_r4(:,:,:,:)
597 real(4), pointer :: cld_ptr_r4(:,:,:,:)
598 real(8), pointer :: pressure(:,:,:)
599 real(8) :: hu, hu_modified
600 real(8) :: cld, cld_modified
601 real(8) :: minValueCld, maxValueCld
602 integer :: lon1, lon2, lat1, lat2, varNameIndex
603 integer :: lonIndex, latIndex, levIndex, stepIndex, varLevIndex, memberIndex
604 integer :: numMember, numStep, numLev, numLev_rttov
605 integer :: fnom, fclos, ierr, nulfile
606 character(len=256) :: fileName
607 character(len=4) :: varName
608 logical, save :: firstTime=.true.
609 logical :: applyLimitToAllVarname
610 logical :: applyLimitToCloud
611
612 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_ens: STARTING'
613
614 if (present(varName_opt)) then
615 applyLimitToAllVarname = .false.
616 varName = varName_opt
617 else
618 applyLimitToAllVarname = .true.
619 varName = 'XXXX'
620 end if
621
622 ! for cloud, limits are applied to ALL cloud variables.
623 if (present(applyLimitToCloud_opt)) then
624 applyLimitToCloud = applyLimitToCloud_opt
625 if (.not. applyLimitToCloud) call utl_abort('qlim_rttovLimit_ens: remove applyLimitToCloud_opt argument')
626 else
627 if (vnl_isCloudVar(varName)) then
628 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_ens: limits are applied to ALL cloud variables'
629 applyLimitToCloud = .true.
630 else
631 applyLimitToCloud = .false.
632 end if
633 end if
634
635 if (applyLimitToCloud) applyLimitToAllVarname = .false.
636
637 if ((applyLimitToAllVarname .or. trim(varName) == 'HU') .and. &
638 ens_varExist(ensemble,'HU')) then
639
640 if ( mmpi_myid == 0 ) write(*,*) 'qlim_rttovLimit_ens: applying limits to HU.'
641
642 if (ens_getDataKind(ensemble) == 8) then
643 call utl_abort('qlim_rttovLimit_ens: Not compatible with dataKind = 8')
644 end if
645
646 ! Read in RTTOV humidity limits
647 fileName = "rttov_h2o_limits.dat"
648 nulfile = 0
649 ierr = fnom(nulfile, fileName, "FMT+OLD+R/O", 0)
650 if( ierr /= 0 ) then
651 if ( mmpi_myid == 0 ) write(*,*) 'fileName = ', fileName
652 call utl_abort('qlim_rttovLimit_ens: error opening the humidity limits file')
653 end if
654
655 read(nulfile,*) numLev_rttov
656 if ( mmpi_myid == 0 .and. firstTime ) write(*,*) 'qlim_rttovLimit_ens: rttov number of levels = ', numLev_rttov
657 allocate(press_rttov(numLev_rttov))
658 allocate(qmin_rttov(numLev_rttov))
659 allocate(qmax_rttov(numLev_rttov))
660 do levIndex = 1, numLev_rttov
661 read(nulfile,*) press_rttov(levIndex), qmax_rttov(levIndex), qmin_rttov(levIndex)
662 end do
663 ierr = fclos(nulfile)
664 press_rttov(:) = press_rttov(:) * mpc_pa_per_mbar_r8
665 qmin_rttov(:) = qmin_rttov(:) / mixratio_to_ppmv
666 qmax_rttov(:) = qmax_rttov(:) / mixratio_to_ppmv
667
668 if (firstTime) then
669 write(*,*) ' '
670 do levIndex = 1, numLev_rttov
671 if ( mmpi_myid == 0 ) write(*,fmt='(" qlim_rttovLimit_ens: LEVEL = ",I4,", PRES = ",F9.0,", HUMIN = ",E10.2,", HUMAX = ",E10.2)') &
672 levIndex, press_rttov(levIndex), qmin_rttov(levIndex), qmax_rttov(levIndex)
673 end do
674 firstTime = .false.
675 end if
676
677 vco_ptr => ens_getVco(ensemble)
678 numLev = ens_getNumLev(ensemble,'TH')
679 numMember = ens_getNumMembers(ensemble)
680 numStep = ens_getNumStep(ensemble)
681 call ens_getLatLonBounds(ensemble, lon1, lon2, lat1, lat2)
682
683 allocate( psfc(numMember,numStep) )
684 if (vco_ptr%vcode == 5100) allocate( psfcLS(numMember,numStep) )
685 allocate( qmin3D_rttov(numMember,numStep,numLev) )
686 allocate( qmax3D_rttov(numMember,numStep,numLev) )
687
688 do latIndex = lat1, lat2
689 do lonIndex = lon1, lon2
690
691 varLevIndex = ens_getKFromLevVarName(ensemble, 1, 'P0')
692 psfc_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
693 psfc(:,:) = real(psfc_ptr_r4(:,:,lonIndex,latIndex),8)
694 if (vco_ptr%vcode == 5100) then
695 varLevIndex = ens_getKFromLevVarName(ensemble, 1, 'P0LS')
696 psfcLS_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
697 psfcLS(:,:) = real(psfcLS_ptr_r4(:,:,lonIndex,latIndex),8)
698 call czp_fetch3DLevels(vco_ptr, psfc, sfcFldLS_opt=psfcLS, fldT_opt=pressure)
699 else
700 call czp_fetch3DLevels(vco_ptr, psfc, fldT_opt=pressure)
701 end if
702
703 ! Interpolate RTTOV limits onto model levels
704 call qlim_lintv_minmax(press_rttov, qmin_rttov, qmax_rttov, numLev_rttov, &
705 numMember, numStep, numLev, pressure, &
706 qmin3D_rttov, qmax3D_rttov)
707
708 do levIndex = 1, numLev
709
710 varLevIndex = ens_getKFromLevVarName(ensemble, levIndex, 'HU')
711 hu_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
712
713 !$OMP PARALLEL DO PRIVATE (stepIndex, memberIndex, hu, hu_modified)
714 do stepIndex = 1, numStep
715 do memberIndex = 1, numMember
716
717 hu = real(hu_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex),8)
718
719 ! limit the humidity according to the rttov limits
720 hu_modified = max(hu, qmin3D_rttov(memberIndex, stepIndex, levIndex) )
721 hu_modified = min(hu_modified, qmax3D_rttov(memberIndex, stepIndex, levIndex) )
722 hu_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = real(hu_modified,4)
723
724 end do ! memberIndex
725 end do ! stepIndex
726 !$OMP END PARALLEL DO
727 end do ! levIndex
728
729 deallocate(pressure)
730
731 end do ! lonIndex
732 end do ! latIndex
733
734 deallocate( psfc )
735 if (allocated(psfcLS)) deallocate( psfcLS )
736 deallocate( qmin3D_rttov )
737 deallocate( qmax3D_rttov )
738
739 deallocate(qmax_rttov)
740 deallocate(qmin_rttov)
741 deallocate(press_rttov)
742 end if
743
744 ! apply limits to ALL available cloud variables
745 if ((applyLimitToAllVarname .or. applyLimitToCloud) .and. cloudExistInEnsemble(ensemble)) then
746
747 vco_ptr => ens_getVco(ensemble)
748 numLev = ens_getNumLev(ensemble,'TH')
749 numMember = ens_getNumMembers(ensemble)
750 numStep = ens_getNumStep(ensemble)
751 call ens_getLatLonBounds(ensemble, lon1, lon2, lat1, lat2)
752
753 do varNameIndex = 1, vnl_numvarmaxCloud
754 if (.not. ens_varExist(ensemble, vnl_varNameListCloud(varNameIndex))) cycle
755
756 if (mmpi_myid == 0) write(*,*) 'qlim_rttovLimit_ens: applying limits to ', &
757 vnl_varNameListCloud(varNameIndex)
758
759 minValueCld = qlim_getMinValueCloud(vnl_varNameListCloud(varNameIndex))
760 maxValueCld = qlim_getMaxValueCloud(vnl_varNameListCloud(varNameIndex))
761
762 do latIndex = lat1, lat2
763 do lonIndex = lon1, lon2
764
765 do levIndex = 1, numLev
766 varLevIndex = ens_getKFromLevVarName(ensemble, levIndex, vnl_varNameListCloud(varNameIndex))
767 cld_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
768
769 !$OMP PARALLEL DO PRIVATE (stepIndex, memberIndex, cld, cld_modified)
770 do stepIndex = 1, numStep
771 do memberIndex = 1, numMember
772
773 cld = real(cld_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex),8)
774
775 cld_modified = max(cld,minValueCld)
776 cld_modified = min(cld_modified,maxValueCld)
777 cld_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = real(cld_modified,4)
778
779 end do ! memberIndex
780 end do ! stepIndex
781 !$OMP END PARALLEL DO
782
783 end do ! levIndex
784
785 end do ! lonIndex
786 end do ! latIndex
787 end do ! varNameIndex
788
789 end if
790
791 end subroutine qlim_rttovLimit_ens
792
793 !--------------------------------------------------------------------------
794 ! qlim_lintv_minmax
795 !--------------------------------------------------------------------------
796 subroutine qlim_lintv_minmax(press_src, qmin_src, qmax_src, numLev_src, &
797 ni_dest, nj_dest, numLev_dest, press_dest, &
798 qmin_dest, qmax_dest)
799 !
800 !:Purpose: To perform the vertical interpolation in log of pressure and
801 ! and constant value extrapolation of one-dimensional vectors.
802 !
803 implicit none
804
805 ! Arguments:
806 real(8), intent(in) :: press_src(numLev_src) ! Vertical levels, pressure (source)
807 real(8), intent(in) :: qmin_src(numLev_src) ! Vectors to be interpolated (source)
808 real(8), intent(in) :: qmax_src(numLev_src) ! Vectors to be interpolated (source)
809 integer, intent(in) :: numLev_src ! Number of input levels (source)
810 integer, intent(in) :: ni_dest ! Number of profiles
811 integer, intent(in) :: nj_dest ! Number of profiles
812 integer, intent(in) :: numLev_dest ! Number of output levels (destination)
813 real(8), intent(in) :: press_dest(:,:,:) ! Vertical levels, pressure (destination)
814 real(8), intent(out) :: qmin_dest(:,:,:) ! Interpolated profiles (destination)
815 real(8), intent(out) :: qmax_dest(:,:,:) ! Interpolated profiles (destination)
816
817 ! Locals:
818 integer :: ji, jk, jo, ii, jj, ik, iorder
819 real(8) :: zpo(numLev_dest)
820 integer :: il(numLev_dest)
821 real(8) :: zpi(0:numLev_src+1)
822 real(8) :: zqmin_src(0:numLev_src+1), zqmax_src(0:numLev_src+1)
823 real(8) :: zw1, zw2, zp, xi, zrt, zp1, zp2
824
825 zpi(0)=200000.d0
826 zpi(numLev_src+1)=200000.d0
827
828 ! Determine if input pressure levels are in ascending or
829 ! descending order.
830 if ( press_src(1) < press_src(numLev_src) ) then
831 iorder = 1
832 else
833 iorder = -1
834 end if
835
836 ! Source levels
837 do jk = 1, numLev_src
838 zpi(jk) = press_src(jk)
839 zqmin_src(jk) = qmin_src(jk)
840 zqmax_src(jk) = qmax_src(jk)
841 enddo
842 zqmin_src(0) = qmin_src(1)
843 zqmin_src(numLev_src+1) = qmin_src(numLev_src)
844 zqmax_src(0) = qmax_src(1)
845 zqmax_src(numLev_src+1) = qmax_src(numLev_src)
846
847 do jj = 1, nj_dest
848 do ii = 1, ni_dest
849 zpo(:) = press_dest(ii,jj,:)
850
851 ! Interpolate in log of pressure or extrapolate with constant value
852 ! for each destination pressure level
853
854 ! Find the adjacent level below
855 il(:) = 0
856 do ji = 1, numLev_src
857 do jo = 1, numLev_dest
858 zrt = zpo(jo)
859 zp = zpi(ji)
860 xi = sign(1.0d0,iorder*(zrt-zp))
861 il(jo) = il(jo) + max(0.0d0,xi)
862 end do
863 end do
864
865 ! Interpolation/extrapolation
866 do jo = 1, numLev_dest
867 ik = il(jo)
868 zp = zpo(jo)
869 zp1 = zpi(ik)
870 zp2 = zpi(ik+1)
871 zw1 = log(zp/zp2)/log(zp1/zp2)
872 zw2 = 1.d0 - zw1
873 qmin_dest(ii,jj,jo) = zw1*zqmin_src(ik) + zw2*zqmin_src(ik+1)
874 qmax_dest(ii,jj,jo) = zw1*zqmax_src(ik) + zw2*zqmax_src(ik+1)
875 end do
876
877 end do ! ii
878 end do ! jj
879
880 end subroutine qlim_lintv_minmax
881
882 !--------------------------------------------------------------------------
883 ! qlim_setMin_ens
884 !--------------------------------------------------------------------------
885 subroutine qlim_setMin_ens(ensemble,huMinValue)
886 !
887 !:Purpose: To impose lower limit on humidity variable of an ensemble
888 !
889 implicit none
890
891 ! Arguments:
892 type(struct_ens), intent(inout) :: ensemble
893 real(8), intent(in) :: huMinValue
894
895 ! Locals:
896 real(4), pointer :: hu_ptr_r4(:,:,:,:)
897 real(4) :: hu, hu_modified
898 integer :: lon1, lon2, lat1, lat2, numLev
899 integer :: lonIndex, latIndex, levIndex, stepIndex, memberIndex, varLevIndex
900
901 if (mmpi_myid == 0) write(*,*) 'qlim_setMin_ens: STARTING'
902
903 if (ens_getDataKind(ensemble) == 8) then
904 call utl_abort('qlim_setMin_ens: Not compatible with dataKind = 8')
905 end if
906
907 if( .not. ens_varExist(ensemble,'HU') ) then
908 if( mmpi_myid == 0 ) write(*,*) 'qlim_setMin_ens: ensemble does not ' // &
909 'contain humidity ... doing nothing'
910 return
911 end if
912
913 numLev = ens_getNumLev(ensemble,'TH')
914 call ens_getLatLonBounds(ensemble, lon1, lon2, lat1, lat2)
915
916 do latIndex = lat1, lat2
917 do lonIndex = lon1, lon2
918 do levIndex = 1, numLev
919 varLevIndex = ens_getKFromLevVarName(ensemble, levIndex, 'HU')
920 hu_ptr_r4 => ens_getOneLev_r4(ensemble,varLevIndex)
921
922 !$OMP PARALLEL DO PRIVATE (stepIndex, memberIndex, hu, hu_modified)
923 do stepIndex = 1, ens_getNumStep(ensemble)
924 do memberIndex = 1, ens_getNumMembers(ensemble)
925 hu = hu_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex)
926 hu_modified = max(hu, real(huMinValue,4))
927 hu_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = hu_modified
928 end do ! memberIndex
929 end do ! stepIndex
930 !$OMP END PARALLEL DO
931
932 end do ! levIndex
933 end do ! lonIndex
934 end do ! latIndex
935
936 end subroutine qlim_setMin_ens
937
938 !-----------------------------------------------------------------------
939 ! qlim_getMinValueCloud
940 !----------------------------------------------------------------------
941 function qlim_getMinValueCloud(varName) result(minValue)
942 !
943 !:Purpose: Return the minValue for the hydrometeor.
944 !
945 implicit none
946
947 ! Arguments:
948 character(len=*), intent(in) :: varName
949 ! Result:
950 real(8) :: minValue
951
952 ! readNameList runs one time during program execution
953 call readNameList
954
955 select case (trim(varName))
956 case ('LWCR')
957 minValue = qlim_minValueLWCR
958 case ('IWCR')
959 minValue = qlim_minValueIWCR
960 case ('RF')
961 minValue = qlim_minValueRF
962 case ('SF')
963 minValue = qlim_minValueSF
964 case ('CLDR')
965 minValue = qlim_minValueCLDR
966 case default
967 write(*,*)
968 write(*,*) 'ERROR unknown varName: ', trim(varName)
969 call utl_abort('qlim_getMinValueCloud')
970 end select
971
972 end function qlim_getMinValueCloud
973
974 !-----------------------------------------------------------------------
975 ! qlim_getMaxValueCloud
976 !----------------------------------------------------------------------
977 function qlim_getMaxValueCloud(varName) result(maxValue)
978 !
979 !:Purpose: Return the maxValue for the hydrometeor.
980 !
981 implicit none
982
983 ! Arguments:
984 character(len=*), intent(in) :: varName
985 ! Result:
986 real(8) :: maxValue
987
988 ! readNameList runs one time during program execution
989 call readNameList
990
991 select case (trim(varName))
992 case ('LWCR')
993 maxValue = qlim_maxValueLWCR
994 case ('IWCR')
995 maxValue = qlim_maxValueIWCR
996 case ('RF')
997 maxValue = qlim_maxValueRF
998 case ('SF')
999 maxValue = qlim_maxValueSF
1000 case ('CLDR')
1001 maxValue = qlim_maxValueCLDR
1002 case default
1003 write(*,*)
1004 write(*,*) 'ERROR unknown varName: ', trim(varName)
1005 call utl_abort('qlim_getMaxValueCloud')
1006 end select
1007
1008 end function qlim_getMaxValueCloud
1009
1010 !-----------------------------------------------------------------------
1011 ! cloudExistInEnsemble
1012 !----------------------------------------------------------------------
1013 function cloudExistInEnsemble(ensemble) result(cloudExist)
1014 !
1015 !:Purpose: determine if any cloud variable exists in the ensemble.
1016 !
1017 implicit none
1018
1019 ! Arguments:
1020 type(struct_ens), intent(in) :: ensemble
1021 ! Result:
1022 logical :: cloudExist
1023
1024 ! Locals:
1025 integer :: varNameIndex
1026
1027 cloudExist = .false.
1028 do varNameIndex = 1, vnl_numvarmaxCloud
1029 if (ens_varExist(ensemble, vnl_varNameListCloud(varNameIndex))) then
1030 cloudExist = .true.
1031 return
1032 end if
1033 end do
1034
1035 end function cloudExistInEnsemble
1036
1037 !-----------------------------------------------------------------------
1038 ! cloudExistInStateVector
1039 !----------------------------------------------------------------------
1040 function cloudExistInStateVector(stateVector) result(cloudExist)
1041 !
1042 !:Purpose: determine if any cloud variable exists in the stateVector.
1043 !
1044 implicit none
1045
1046 ! Arguments:
1047 type(struct_gsv), intent(in) :: stateVector
1048 ! Result:
1049 logical :: cloudExist
1050
1051 ! Locals:
1052 integer :: varNameIndex
1053
1054 cloudExist = .false.
1055 do varNameIndex = 1, vnl_numvarmaxCloud
1056 if (gsv_varExist(stateVector, vnl_varNameListCloud(varNameIndex))) then
1057 cloudExist = .true.
1058 return
1059 end if
1060 end do
1061
1062 end function cloudExistInStateVector
1063
1064end module humidityLimits_mod
