UPP (develop)
Loading...
Searching...
No Matches
CALRAD_WCLOUD_newcrtm.f
Go to the documentation of this file.
1
25!---------------------------------------------------------------------------
27!---------------------------------------------------------------------------
28 SUBROUTINE calrad_wcloud
29
30 use vrbls3d, only: o3, pint, pmid, t, q, qqw, qqi, qqr, f_rimef, nlice, nrain, qqs, qqg, &
31 qqnr, qqni, qqnw, cfr, effri, effrl, effrs
32 use vrbls2d, only: czen, ivgtyp, sno, pctsno, ths, vegfrc, si, u10h, v10h, u10,&
33 v10, smstot, hbot, htop, cnvcfr
34 use masks, only: gdlat, gdlon, sm, lmh, sice
35 use soil, only:
36 use gridspec_mod, only: gridtype
37 use cmassi_mod, only: trad_ice
38 use kinds, only: r_kind,r_single,r_double,i_kind
39 use crtm_module, only: crtm_atmosphere_type,crtm_surface_type,crtm_geometry_type, &
40 crtm_surface_create,o3_id,co2_id,crtm_forward,mass_mixing_ratio_units, &
41 crtm_atmosphere_create, &
42 crtm_options_type,crtm_destroy,crtm_init,specific_amount_units, &
43 success,crtm_options_destroy,crtm_options_create, crtm_options_associated
44
45 use crtm_rtsolution_define, only: crtm_rtsolution_type, crtm_rtsolution_create, &
46 crtm_rtsolution_destroy, crtm_rtsolution_associated
47 use crtm_spccoeff, only: sc
48 use crtm_atmosphere_define, only:h2o_id,crtm_atmosphere_associated, &
49 crtm_atmosphere_destroy,volume_mixing_ratio_units,crtm_atmosphere_zero
50 use crtm_surface_define, only: crtm_surface_destroy, crtm_surface_associated, &
51 crtm_surface_zero
52 use crtm_channelinfo_define, only: crtm_channelinfo_type
53! use crtm_channelinfo_define, only: crtm_channelinfo_type, AntCorr_type
54 use crtm_parameters, only: limit_exp,toa_pressure,max_n_layers,max_sensor_scan_angle
55 use crtm_cloud_define, only: water_cloud,ice_cloud,rain_cloud,snow_cloud,graupel_cloud,hail_cloud
56 use message_handler, only: success,warning, display_message
57
58 use params_mod, only: pi, rtd, p1000, capa, h1000, h1, g, rd, d608, qconv, small
59 use rqstfld_mod, only: iget, id, lvls, iavblfld
60 use ctlblk_mod, only: modelname, ivegsrc, novegtype, imp_physics, lm, spval, icu_physics,&
61 grib, cfld, fld_info, datapd, idat, im, jsta, jend, jm, me, ista, iend
62!
63 implicit none
64
65 ! DECLARE VARIABLES.
66 !
67 ! Mapping land surface type of GFS to CRTM
68 ! Note: index 0 is water, and index 13 is ice. The two indices are not
69 ! used and just assigned to COMPACTED_SOIL.
70 !integer, parameter, dimension(0:13) :: gfs_to_crtm=(/COMPACTED_SOIL, &
71 ! BROADLEAF_FOREST, BROADLEAF_FOREST, BROADLEAF_PINE_FOREST, PINE_FOREST, &
72 ! PINE_FOREST, BROADLEAF_BRUSH, SCRUB, SCRUB, SCRUB_SOIL, TUNDRA, &
73 ! COMPACTED_SOIL, TILLED_SOIL, COMPACTED_SOIL/)
74
75 ! Mapping land surface type of NMM to CRTM
76 ! Note: index 16 is water, and index 24 is ice. The two indices are not
77 ! used and just assigned to COMPACTED_SOIL.
78 ! integer, parameter, dimension(24) :: nmm_to_crtm=(/URBAN_CONCRETE, &
79 ! & COMPACTED_SOIL, IRRIGATED_LOW_VEGETATION, GRASS_SOIL, MEADOW_GRASS, &
80 ! & MEADOW_GRASS, MEADOW_GRASS, SCRUB, GRASS_SCRUB, MEADOW_GRASS, &
81 ! & BROADLEAF_FOREST, PINE_FOREST, BROADLEAF_FOREST, PINE_FOREST, &
82 ! & BROADLEAF_PINE_FOREST, COMPACTED_SOIL, WET_SOIL, WET_SOIL, &
83 ! & IRRIGATED_LOW_VEGETATION, TUNDRA, TUNDRA, TUNDRA, TUNDRA, &
84 ! & COMPACTED_SOIL/)
85
86 ! For land, the land types
87 !INTEGER, PARAMETER :: N_VALID_LAND_TYPES = 20
88 INTEGER, PARAMETER :: INVALID_LAND = 0
89 INTEGER, PARAMETER :: COMPACTED_SOIL = 1
90 INTEGER, PARAMETER :: TILLED_SOIL = 2
91 INTEGER, PARAMETER :: SAND = 3
92 INTEGER, PARAMETER :: ROCK = 4
93 INTEGER, PARAMETER :: IRRIGATED_LOW_VEGETATION = 5
94 INTEGER, PARAMETER :: MEADOW_GRASS = 6
95 INTEGER, PARAMETER :: SCRUB = 7
96 INTEGER, PARAMETER :: BROADLEAF_FOREST = 8
97 INTEGER, PARAMETER :: PINE_FOREST = 9
98 INTEGER, PARAMETER :: TUNDRA = 10
99 INTEGER, PARAMETER :: GRASS_SOIL = 11
100 INTEGER, PARAMETER :: BROADLEAF_PINE_FOREST = 12
101 INTEGER, PARAMETER :: GRASS_SCRUB = 13
102 INTEGER, PARAMETER :: SOIL_GRASS_SCRUB = 14
103 INTEGER, PARAMETER :: URBAN_CONCRETE = 15
104 INTEGER, PARAMETER :: PINE_BRUSH = 16
105 INTEGER, PARAMETER :: BROADLEAF_BRUSH = 17
106 INTEGER, PARAMETER :: WET_SOIL = 18
107 INTEGER, PARAMETER :: SCRUB_SOIL = 19
108 INTEGER, PARAMETER :: BROADLEAF70_PINE30 = 20
109
110
111 integer, allocatable:: model_to_crtm(:)
112 integer, parameter:: ndat=100
113 ! CRTM structure variable declarations.
114 integer,parameter:: n_absorbers = 2
115 ! integer,parameter:: n_clouds = 4
116 integer,parameter:: n_aerosols = 0
117 ! Add your sensors here
118 integer(i_kind),parameter:: n_sensors=23
119 character(len=20),parameter,dimension(1:n_sensors):: sensorlist= &
120 (/'imgr_g15 ', &
121 'imgr_g13 ', &
122 'imgr_g12 ', &
123 'imgr_g11 ', &
124 'amsre_aqua ', &
125 'tmi_trmm ', &
126 'ssmi_f13 ', &
127 'ssmi_f14 ', &
128 'ssmi_f15 ', &
129 'ssmis_f16 ', &
130 'ssmis_f17 ', &
131 'ssmis_f18 ', &
132 'ssmis_f19 ', &
133 'ssmis_f20 ', &
134 'seviri_m10 ', &
135 'imgr_mt2 ', &
136 'imgr_mt1r ', &
137 'imgr_insat3d ', &
138 'abi_gr ', &
139 'abi_g16 ', &
140 'abi_g17 ', &
141 'abi_g18 ', &
142 'ahi_himawari8 '/)
143 character(len=13),parameter,dimension(1:n_sensors):: obslist= &
144 (/'goes_img ', &
145 'goes_img ', &
146 'goes_img ', &
147 'goes_img ', &
148 'amsre ', &
149 'tmi ', &
150 'ssmi ', &
151 'ssmi ', &
152 'ssmi ', &
153 'ssmis ', &
154 'ssmis ', &
155 'ssmis ', &
156 'ssmis ', &
157 'ssmis ', &
158 'seviri ', &
159 'imgr_mt2 ', &
160 'imgr_mt1r ', &
161 'imgr_insat3d ', &
162 'abi ', &
163 'abi ', &
164 'abi ', &
165 'abi ', &
166 'ahi_himawari8'/)
167 character(len=20),dimension(1:n_sensors):: sensorlist_local
168!
169 integer(i_kind) sensorindex
170 integer(i_kind) lunin,nobs,nchanl,nreal
171 integer(i_kind) error_status,itype
172 integer(i_kind) err1,err2,err3,err4
173 integer(i_kind) i,j,k,msig
174 integer(i_kind) lcbot,lctop !bsf
175 integer jdn,ichan,ixchan,igot
176 integer isat
177
178! Wm Lewis: added
179 real :: EFFR
180
181 real(r_kind),parameter:: r100=100.0_r_kind
182 real,parameter:: ozsmall = 1.e-10 ! to convert to mass mixing ratio
183 real(r_kind) tsfc
184 real(r_double),dimension(4):: sfcpct
185 real(r_kind) snodepth,vegcover
186 real snoeqv
187 real snofrac
188 real(r_kind),dimension(ista:iend,jsta:jend):: tb1,tb2,tb3,tb4
189 real(r_kind),allocatable :: tb(:,:,:)
190 real,dimension(im,jm):: grid1
191 real sun_zenith,sun_azimuth, dpovg, sun_zenith_rad
192 real sat_zenith
193 real q_conv !bsf
194 real,parameter:: constoz = 604229.0_r_kind
195 real sublat,sublon
196 real RHO,RHOX
197 character(13)::obstype
198 character(20)::isis
199 character(20)::isis_local
200
201 logical hirs2,msu,goessndr,hirs3,hirs4,hirs,amsua,amsub,airs,hsb &
202 ,goes_img,abi,seviri, mhs,insat3d
203 logical avhrr,avhrr_navy,lextra,ssu
204 logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,change
205 logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img
206 logical sea,mixed,land,ice,snow,toss
207 logical micrim,microwave
208 logical post_abig16, post_abig17, post_abig18, post_abigr ! if true, user requested at least one abi channel
209 logical fix_abig16, fix_abig17, fix_abig18 ! if true, abi_g16, abi_g17 fix files are available
210 logical post_ahi8 ! if true, user requested at least on ahi channel (himawari8)
211 logical post_ssmis17 ! if true, user requested at least on ssmis_f17 channel
212 ! logical,dimension(nobs):: luse
213 logical, parameter :: debugprint = .false.
214 type(crtm_atmosphere_type),dimension(1):: atmosphere
215 type(crtm_surface_type),dimension(1) :: surface
216 type(crtm_geometry_type),dimension(1) :: geometryinfo
217 type(crtm_options_type),dimension(1) :: options
218
219 type(crtm_rtsolution_type),allocatable,dimension(:,:):: rtsolution
220 type(crtm_channelinfo_type),allocatable,dimension(:) :: channelinfo
221!
222 integer ii,jj,n_clouds,n,nc
223 integer,external :: iw3jdn
224 !
225
226 !*****************************************************************************
227 ! This code and sensorlist_local, isis_local can be modified/removed when the
228 ! linked CRTM version is updated with fix files abi_g16 & abi_g17
229 fix_abig16 = .false.
230 fix_abig17 = .false.
231 fix_abig18 = .false.
232 do n=1, n_sensors
233 sensorlist_local(n) = sensorlist(n)
234 if (sensorlist(n) == 'abi_g16') then ! check if fix file is available
235 inquire(file='abi_g16.SpcCoeff.bin',exist=fix_abig16)
236 if (.not.fix_abig16) sensorlist_local(n) = 'abi_gr '
237 endif
238 if (sensorlist(n) == 'abi_g17') then
239 inquire(file='abi_g17.SpcCoeff.bin',exist=fix_abig17)
240 if (.not.fix_abig17) sensorlist_local(n) = 'abi_gr '
241 endif
242 if (sensorlist(n) == 'abi_g18') then
243 inquire(file='abi_g18.SpcCoeff.bin',exist=fix_abig18)
244 if (.not.fix_abig18) sensorlist_local(n) = 'abi_gr '
245 endif
246 enddo
247
248
249 ! Mapping land surface type of NMM to CRTM
250 !if(MODELNAME == 'NMM' .OR. MODELNAME == 'NCAR' .OR. MODELNAME == 'RAPR')then
251 if(ivegsrc==1)then !IGBP veg type
252 allocate(model_to_crtm(novegtype) )
253 model_to_crtm=(/pine_forest, broadleaf_forest, pine_forest, &
254 broadleaf_forest,broadleaf_pine_forest, scrub, scrub_soil, &
255 broadleaf_brush,broadleaf_brush, scrub, broadleaf_brush, &
256 tilled_soil, urban_concrete,tilled_soil, urban_concrete, &
257 compacted_soil, broadleaf_brush, tundra,tundra, tundra/)
258 else if(ivegsrc==0)then ! USGS veg type
259 allocate(model_to_crtm(novegtype) )
260 model_to_crtm=(/urban_concrete, &
261 compacted_soil, irrigated_low_vegetation, grass_soil, meadow_grass, &
262 meadow_grass, meadow_grass, scrub, grass_scrub, meadow_grass, &
263 broadleaf_forest, pine_forest, broadleaf_forest, pine_forest, &
264 broadleaf_pine_forest, compacted_soil, wet_soil, wet_soil, &
265 irrigated_low_vegetation, tundra, tundra, tundra, tundra, &
266 compacted_soil/)
267 else if(ivegsrc==2)then ! old GFS veg type
268 allocate(model_to_crtm(0:novegtype) )
269 model_to_crtm=(/compacted_soil, &
270 broadleaf_forest, broadleaf_forest, broadleaf_pine_forest, &
271 pine_forest, pine_forest, broadleaf_brush, scrub, scrub, scrub_soil, &
272 tundra, compacted_soil, tilled_soil, compacted_soil/)
273 else
274 print*,'novegtype=',novegtype
275 print*,'model veg type not supported by post in calling crtm '
276 print*,'skipping generation of simulated radiance'
277 return
278 end if
279 !end if
280
281 !10 channels, easier to set a logical
282 post_abig16=.false.
283 do n = 927, 927+9 ! 927 set in RQSTFLD.f
284 if (iget(n) > 0) post_abig16=.true.
285 enddo
286 post_abig17=.false.
287 do n = 937, 937+9 ! 937 set in RQSTFLD.f
288 if (iget(n) > 0) post_abig17=.true.
289 enddo
290 post_abig18=.false.
291 do n = 531, 531+9 ! 531 set in RQSTFLD.f
292 if (iget(n) > 0) post_abig18=.true.
293 enddo
294 post_abigr=.false.
295 do n = 958, 958+9 ! 958 set in RQSTFLD.f
296 if (iget(n) > 0) post_abigr=.true.
297 enddo
298 post_ahi8=.false.
299 do n = 969, 969+9 ! 969 set in RQSTFLD.f
300 if (iget(n) > 0) post_ahi8=.true.
301 enddo
302 post_ssmis17=.false.
303 do n = 825, 825+3 ! 825 set in RQSTFLD.f
304 if (iget(n) > 0) post_ssmis17=.true.
305 enddo
306
307 ! DO NOT FORGET TO ADD YOUR NEW IGET HERE (IF YOU'VE ADDED ONE)
308 ! START SUBROUTINE CALRAD.
309 ifactive: if (iget(327) > 0 .or. iget(328) > 0 .or. iget(329) > 0 &
310 .or. iget(330) > 0 .or. iget(446) > 0 .or. iget(447) > 0 &
311 .or. iget(448) > 0 .or. iget(449) > 0 .or. iget(456) > 0 &
312 .or. iget(457) > 0 .or. iget(458) > 0 .or. iget(459) > 0 &
313 .or. iget(460) > 0 .or. iget(461) > 0 .or. iget(462) > 0 &
314 .or. iget(463) > 0 .or. iget(483) > 0 .or. iget(484) > 0 &
315 .or. iget(485) > 0 .or. iget(486) > 0 .or. iget(488) > 0 &
316 .or. iget(489) > 0 .or. iget(490) > 0 .or. iget(491) > 0 &
317 .or. iget(492) > 0 .or. iget(493) > 0 .or. iget(494) > 0 &
318 .or. iget(495) > 0 .or. iget(496) > 0 .or. iget(497) > 0 &
319 .or. iget(498) > 0 .or. iget(499) > 0 .or. iget(800) > 0 &
320 .or. iget(801) > 0 .or. iget(802) > 0 .or. iget(803) > 0 &
321 .or. iget(804) > 0 .or. iget(805) > 0 .or. iget(806) > 0 &
322 .or. iget(807) > 0 .or. iget(809) > 0 &
323 .or. iget(810) > 0 .or. iget(811) > 0 .or. iget(812) > 0 &
324 .or. iget(813) > 0 .or. iget(814) > 0 .or. iget(815) > 0 &
325 .or. iget(816) > 0 .or. iget(817) > 0 .or. iget(818) > 0 &
326 .or. iget(819) > 0 .or. iget(820) > 0 .or. iget(821) > 0 &
327 .or. iget(822) > 0 .or. iget(823) > 0 .or. iget(824) > 0 &
328 .or. iget(825) > 0 .or. iget(826) > 0 .or. iget(827) > 0 &
329 .or. iget(828) > 0 .or. iget(829) > 0 .or. iget(830) > 0 &
330 .or. iget(831) > 0 .or. iget(832) > 0 .or. iget(833) > 0 &
331 .or. iget(834) > 0 .or. iget(835) > 0 .or. iget(836) > 0 &
332 .or. iget(837) > 0 .or. iget(838) > 0 .or. iget(839) > 0 &
333 .or. iget(840) > 0 .or. iget(841) > 0 .or. iget(842) > 0 &
334 .or. iget(843) > 0 .or. iget(844) > 0 .or. iget(845) > 0 &
335 .or. iget(846) > 0 .or. iget(847) > 0 .or. iget(848) > 0 &
336 .or. iget(849) > 0 .or. iget(850) > 0 .or. iget(851) > 0 &
337 .or. iget(852) > 0 .or. iget(856) > 0 .or. iget(857) > 0 &
338 .or. iget(860) > 0 .or. iget(861) > 0 &
339 .or. iget(862) > 0 .or. iget(863) > 0 .or. iget(864) > 0 &
340 .or. iget(865) > 0 .or. iget(866) > 0 .or. iget(867) > 0 &
341 .or. iget(868) > 0 .or. iget(869) > 0 .or. iget(870) > 0 &
342 .or. iget(871) > 0 .or. iget(872) > 0 .or. iget(873) > 0 &
343 .or. iget(874) > 0 .or. iget(875) > 0 .or. iget(876) > 0 &
344 .or. iget(877) > 0 .or. iget(878) > 0 .or. iget(879) > 0 &
345 .or. iget(880) > 0 .or. iget(881) > 0 .or. iget(882) > 0 &
346 .or. post_ahi8 .or. post_ssmis17 &
347 .or. post_abig16 .or. post_abig17 .or. post_abig18 &
348 .or. post_abigr ) then
349
350 ! specify numbers of cloud species
351 ! thompson==8, ferrier==5,95, wsm6==6, lin==2
352 if(imp_physics==99 .or. imp_physics==98)then ! Zhao Scheme
353 n_clouds=2 ! GFS uses Zhao scheme
354 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
355 n_clouds=6 ! change to 6 cloud types because microwave is sensitive to density
356 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 &
357 .or. imp_physics==28 .or. imp_physics==11)then
358 n_clouds=5
359 else
360 n_clouds=0
361 print*,'Warning: number of cloud species (n_clouds) being set to zero for imp_physics=',imp_physics
362 end if
363
364 ! initialize debug print gridpoint index to middle of tile:
365 ii=im/2
366 jj=(jsta+jend)/2
367
368 ! initialize ozone to zeros for wrf nmm and arw for now
369 if (modelname == 'NMM' .OR. modelname == 'NCAR' .OR. modelname == 'RAPR' &
370 )o3=0.0
371 ! compute solar zenith angle for gfs, arw now computes czen in initpost
372! if (MODELNAME == 'GFS')then
373 jdn=iw3jdn(idat(3),idat(1),idat(2))
374 do j=jsta,jend
375 do i=ista,iend
376 call zensun(jdn,float(idat(4)),gdlat(i,j),gdlon(i,j) &
377 ,pi,sun_zenith,sun_azimuth)
378 sun_zenith_rad=sun_zenith/rtd
379 czen(i,j)=cos(sun_zenith_rad)
380 end do
381 end do
382 if(jj>=jsta .and. jj<=jend.and.debugprint) &
383 print*,'sample GFS zenith angle=',acos(czen(ii,jj))*rtd
384! end if
385 ! initialize crtm. load satellite sensor array.
386 ! the optional arguments process_id and output_process_id limit
387 ! generation of runtime informative output to mpi task
388 ! output_process_id(which here is set to be task 0)
389 if(me==0)print*,'success in CALRAD= ',success
390 allocate( channelinfo(n_sensors))
391
392 error_status = crtm_init(sensorlist_local,channelinfo, &
393 process_id=0,output_process_id=0 )
394 if(me==0)print*, 'channelinfo after init= ',channelinfo(1)%sensor_id, &
395 channelinfo(2)%sensor_id
396 if (error_status /= 0_i_kind) &
397 write(6,*)'ERROR*** crtm_init error_status=',error_status
398
399 ! restrict channel list to those which are selected for simulation
400 ! in the lvls filed of wrf_cntrl.parm(does not currently apply
401 ! to all sensors / channels).
402
403 ! goes-13
404 if(iget(868)>0)then
405 call select_channels_l(channelinfo(2),4,(/ 1,2,3,4 /),lvls(1:4,iget(868)),iget(868))
406 endif
407 ! goes-15
408 if(iget(872)>0)then
409 call select_channels_l(channelinfo(1),4,(/ 1,2,3,4 /),lvls(1:4,iget(872)),iget(872))
410 endif
411 ! goes-16
412 if(post_abig16)then
413 nchanl=0
414 do n = 927, 927+9 ! 927 set in RQSTFLD.f
415 if (iget(n) > 0) then
416 nchanl = nchanl+1
417 endif
418 enddo
419 if (nchanl > 0 .and. nchanl <10) then
420 do n = 927, 927+9 ! 927 set in RQSTFLD.f
421 if (iget(n) == 0) channelinfo(19)%Process_Channel(n-927+1)=.false. ! turn off channel processing
422 enddo
423 endif
424 endif
425 ! goes-17
426 if(post_abig17)then
427 nchanl=0
428 do n = 937, 937+9 ! 937 set in RQSTFLD.f
429 if (iget(n) > 0) then
430 nchanl = nchanl+1
431 endif
432 enddo
433 if (nchanl > 0 .and. nchanl <10) then
434 do n = 937, 937+9 ! 927 set in RQSTFLD.f
435 if (iget(n) == 0) channelinfo(20)%Process_Channel(n-937+1)=.false. ! turn off channel processing
436 enddo
437 endif
438 endif
439 ! goes-18
440 if(post_abig18)then
441 nchanl=0
442 do n = 531, 531+9 ! 531 set in RQSTFLD.f
443 if (iget(n) > 0) then
444 nchanl = nchanl+1
445 endif
446 enddo
447 if (nchanl > 0 .and. nchanl <10) then
448 do n = 531, 531+9 ! 531 set in RQSTFLD.f
449 if (iget(n) == 0) channelinfo(21)%Process_Channel(n-531+1)=.false.
450 enddo
451 endif
452 endif
453 ! goes-r for nadir output
454 if(post_abigr)then
455 nchanl=0
456 do n = 958, 958+9 ! 958 set in RQSTFLD.f
457 if (iget(n) > 0) then
458 nchanl = nchanl+1
459 endif
460 enddo
461 if (nchanl > 0 .and. nchanl <10) then
462 do n = 958, 958+9 ! 958 set in RQSTFLD.f
463 if (iget(n) == 0) channelinfo(20)%Process_Channel(n-958+1)=.false. ! turn off channel processing
464 enddo
465 endif
466 endif
467
468 ! himawari-8 ahi infrared
469 if(post_ahi8)then
470 nchanl=0
471 do n = 969, 969+9 ! 969 set in RQSTFLD.f
472 if (iget(n) > 0) then
473 nchanl = nchanl+1
474 endif
475 enddo
476 if (nchanl > 0 .and. nchanl <10) then
477 do n = 969, 969+9 ! 969 set in RQSTFLD.f
478 if (iget(n) == 0) channelinfo(22)%Process_Channel(n-969+1)=.false. ! turn off channel processing
479 enddo
480 endif
481 endif
482
483 ! ssmis f17(37h, 37v, 91h, 91v)
484 if(post_ssmis17)then
485 nchanl=14
486 do n = 825, 825+3 ! 825 set in RQSTFLD.f
487 if (iget(n) > 0) then
488 nchanl = nchanl+1
489 endif
490 enddo
491 if (nchanl > 14 .and. nchanl < 19) then
492 do n = 825, 825+3 ! 825 set in RQSTFLD.f
493 if (iget(n) == 0) channelinfo(11)%Process_Channel(n-825+15)=.false. ! turn off channel processing
494 enddo
495 endif
496 endif
497
498 ! ssmi, f13-f15(19h,19v,??h,37h,37v,85h,85v)
499 if(iget(800)>0)then
500 call select_channels_l(channelinfo(7),7,(/ 1,2,3,4,5,6,7 /),lvls(1:7,iget(800)),iget(800))
501 endif
502 if(iget(806)>0)then
503 call select_channels_l(channelinfo(8),7,(/ 1,2,3,4,5,6,7 /),lvls(1:7,iget(806)),iget(806))
504 endif
505 if(iget(812)>0)then
506 call select_channels_l(channelinfo(9),7,(/ 1,2,3,4,5,6,7 /),lvls(1:7,iget(812)),iget(812))
507 endif
508 ! ssmis, f16-f20(183h,19h,19v,37h,37v,91h,91v)
509 if(iget(818)>0)then
510 call select_channels_l(channelinfo(10),24,(/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24 /),lvls(1:24,iget(818)),iget(818))
511 endif
512! if(iget(825)>0)then
513! call select_channels_L(channelinfo(11),24,(/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24 /),lvls(1:24,iget(825)),iget(825))
514! endif
515 if(iget(832)>0)then
516 call select_channels_l(channelinfo(12),24,(/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24 /),lvls(1:24,iget(832)),iget(832))
517 endif
518 if(iget(839)>0)then
519 call select_channels_l(channelinfo(13),24,(/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24 /),lvls(1:24,iget(839)),iget(839))
520 endif
521 if(iget(846)>0)then
522 call select_channels_l(channelinfo(14),24,(/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24 /),lvls(1:24,iget(846)),iget(846))
523 endif
524 ! seviri
525 if(iget(876)>0)then
526 call select_channels_l(channelinfo(15),8,(/ 1,2,3,4,5,6,7,8 /),lvls(1:8,iget(876)),iget(876))
527 endif
528 ! mt2
529 if(iget(860)>0)then
530 call select_channels_l(channelinfo(16),4,(/ 1,2,3,4 /),lvls(1:4,iget(860)),iget(860))
531 endif
532 ! mt1r
533 if(iget(864)>0)then
534 call select_channels_l(channelinfo(17),4,(/ 1,2,3,4 /),lvls(1:4,iget(864)),iget(864))
535 endif
536 ! insat 3d(kalpana)
537 if(iget(865)>0)then
538 call select_channels_l(channelinfo(18),4,(/ 1,2,3,4 /),lvls(1:4,iget(865)),iget(865))
539 endif
540
541 ! loop over data types to process
542 sensordo: do isat=1,n_sensors
543
544 if(me==0)print*,'n_sensor,obstype,isis',isat,obslist(isat),sensorlist(isat)
545
546 obstype=obslist(isat)
547 isis=trim(sensorlist(isat))
548
549 sensor_avail: if( &
550 (isis=='imgr_g12' .and. (iget(327) > 0 .or. iget(328) > 0 &
551 .or. iget(329) > 0 .or. iget(330) > 0 .or. iget(456) > 0 &
552 .or. iget(457) > 0 .or. iget(458) > 0 .or. iget(459) > 0 )) .OR. &
553 (isis=='imgr_g11' .and. (iget(446) > 0 .or. iget(447) > 0 &
554 .or. iget(448) > 0 .or. iget(449) > 0 .or. iget(460) > 0 &
555 .or. iget(461) > 0 .or. iget(462) > 0 .or. iget(463) > 0)) .OR. &
556 (isis=='amsre_aqua' .and. (iget(483) > 0 .or. iget(484) > 0 &
557 .or. iget(485) > 0 .or. iget(486) > 0)) .OR. &
558 (isis=='tmi_trmm' .and. (iget(488) > 0 .or. iget(489) > 0 &
559 .or. iget(490) > 0 .or. iget(491) > 0)) .OR. &
560 (isis=='ssmi_f13' .and. iget(800) > 0 ) .OR. &
561 (isis=='ssmi_f14' .and. iget(806) > 0 ) .OR. &
562 (isis=='ssmi_f15' .and. iget(812) > 0 ) .OR. &
563 (isis=='ssmis_f16' .and. iget(818) > 0) .OR. &
564 (isis=='ssmis_f17' .and. post_ssmis17) .OR. &
565 (isis=='ssmis_f18' .and. iget(832) > 0) .OR. &
566 (isis=='ssmis_f19' .and. iget(839) > 0) .OR. &
567 (isis=='ssmis_f20' .and. iget(846) > 0) .OR. &
568 (isis=='imgr_mt2' .and. iget(860)>0) .OR. &
569 (isis=='imgr_mt1r' .and. iget(864)>0) .OR. &
570 (isis=='imgr_insat3d' .and. iget(865)>0) .OR. &
571 (isis=='imgr_g13' .and. iget(868)>0) .OR. &
572 (isis=='imgr_g15' .and. iget(872)>0) .OR. &
573 (isis=='abi_g16' .and. post_abig16) .OR. &
574 (isis=='abi_g17' .and. post_abig17) .OR. &
575 (isis=='abi_g18' .and. post_abig18) .OR. &
576 (isis=='abi_gr' .and. post_abigr) .OR. &
577 (isis=='seviri_m10' .and. iget(876)>0) .OR. &
578 (isis=='ahi_himawari8' .and. post_ahi8) )then
579 if(me==0)print*,'obstype, isis= ',obstype,isis
580 ! isis='amsua_n15'
581
582 ! Initialize logical flags for satellite platform
583
584 hirs2 = obstype == 'hirs2'
585 hirs3 = obstype == 'hirs3'
586 hirs4 = obstype == 'hirs4'
587 hirs = hirs2 .or. hirs3 .or. hirs4
588 msu = obstype == 'msu'
589 ssu = obstype == 'ssu'
590 goessndr = obstype == 'sndr' .or. obstype == 'sndrd1' .or. &
591 obstype == 'sndrd2'.or. obstype == 'sndrd3' .or. &
592 obstype == 'sndrd4'
593 amsua = obstype == 'amsua'
594 amsub = obstype == 'amsub'
595 mhs = obstype == 'mhs'
596 airs = obstype == 'airs'
597 hsb = obstype == 'hsb'
598 goes_img = obstype == 'goes_img'
599 abi = obstype == 'abi'
600 seviri = obstype == 'seviri'
601 insat3d = obstype == 'imgr_insat3d'
602 avhrr = obstype == 'avhrr'
603 avhrr_navy = obstype == 'avhrr_navy'
604 ssmi = obstype == 'ssmi'
605 amsre_low = obstype == 'amsre_low'
606 amsre_mid = obstype == 'amsre_mid'
607 amsre_hig = obstype == 'amsre_hig'
608 amsre = amsre_low .or. amsre_mid .or. amsre_hig
609 ssmis = obstype == 'ssmis'
610 ssmis_las = obstype == 'ssmis_las'
611 ssmis_uas = obstype == 'ssmis_uas'
612 ssmis_img = obstype == 'ssmis_img'
613 ssmis_env = obstype == 'ssmis_env'
614
615 ssmis=ssmis_las.or.ssmis_uas.or.ssmis_img.or.ssmis_env.or.ssmis
616
617 micrim=ssmi .or. ssmis .or. amsre ! only used for MW-imager-QC and id_qc(ch)
618
619 microwave=amsua .or. amsub .or. mhs .or. msu .or. hsb .or. micrim
620 ! check sensor list
621 sensorindex = 0
622 sensor_search: do j = 1, n_sensors
623 isis_local = isis ! allows abi_g16 & abi_g17 output using abi_gr fix files
624 if (isis=='abi_g16' .and. .not.fix_abig16) then
625 isis_local='abi_gr '
626 endif
627 if (isis=='abi_g17' .and. .not.fix_abig17) then
628 isis_local='abi_gr '
629 endif
630 if (isis=='abi_g18' .and. .not.fix_abig18) then
631 isis_local='abi_gr '
632 endif
633 if (channelinfo(j)%sensor_id == isis_local ) then
634 sensorindex = j
635 exit sensor_search
636 endif
637 end do sensor_search
638 if (sensorindex == 0 ) then
639 write(6,*)'SETUPRAD: ***WARNING*** problem with sensorindex=',isis,&
640 ' --> CAN NOT PROCESS isis=',isis,' TERMINATE PROGRAM EXECUTION'
641 stop 19
642 endif
643
644! set Satellite IDs for F19 and F20 to valid values since CRTM will not
645! simulate an instrument w/o a WMO ID:
646 if(isis=='ssmis_f19')channelinfo(sensorindex)%WMO_Satellite_Id=287
647 if(isis=='ssmis_f20')channelinfo(sensorindex)%WMO_Satellite_Id=289
648! quiet verbose output warning messages
649 if(isis=='abi_g16')channelinfo(sensorindex)%WMO_Satellite_Id=270
650 if(isis=='abi_g16')channelinfo(sensorindex)%WMO_Sensor_Id=617
651 if(isis=='abi_g17')channelinfo(sensorindex)%WMO_Satellite_Id=271
652 if(isis=='abi_g17')channelinfo(sensorindex)%WMO_Sensor_Id=617
653! assuming sat id for g18 is 272 (continuity w/ g16 and g17)
654 if(isis=='abi_g18')channelinfo(sensorindex)%WMO_Satellite_Id=272
655 if(isis=='abi_g18')channelinfo(sensorindex)%WMO_Sensor_Id=617
656 if(isis=='abi_gr')channelinfo(sensorindex)%WMO_Satellite_Id=270
657 if(isis=='abi_gr')channelinfo(sensorindex)%WMO_Sensor_Id=617
658
659 allocate(rtsolution(channelinfo(sensorindex)%n_channels,1))
660 allocate(tb(ista:iend,jsta:jend,channelinfo(sensorindex)%n_channels))
661 err1=0; err2=0; err3=0; err4=0
662 if(lm > max_n_layers)then
663 write(6,*) 'CALRAD: lm > max_n_layers - '// &
664 'increase crtm max_n_layers ', &
665 lm,max_n_layers
666 stop 2
667 end if
668
669 CALL crtm_atmosphere_create(atmosphere(1),lm,n_absorbers,n_clouds &
670 ,n_aerosols)
671 CALL crtm_surface_create(surface(1),channelinfo(sensorindex)%n_channels)
672 CALL crtm_rtsolution_create(rtsolution,lm)
673 if (.NOT.(crtm_atmosphere_associated(atmosphere(1)))) &
674 write(6,*)' ***ERROR** creating atmosphere.'
675 if (.NOT.(crtm_surface_associated(surface(1)))) &
676 write(6,*)' ***ERROR** creating surface.'
677 if (.NOT.(any(crtm_rtsolution_associated(rtsolution)))) &
678 write(6,*)' ***ERROR** creating rtsolution.'
679
680 atmosphere(1)%n_layers = lm
681 ! atmosphere(1)%level_temperature_input = 0
682 atmosphere(1)%absorber_id(1) = h2o_id
683 atmosphere(1)%absorber_id(2) = o3_id
684 atmosphere(1)%absorber_units(1) = mass_mixing_ratio_units
685 atmosphere(1)%absorber_units(2) = volume_mixing_ratio_units
686 ! atmosphere(1)%absorber_units(2) = MASS_MIXING_RATIO_UNITS ! Use mass mixing ratio
687 atmosphere(1)%level_pressure(0) = toa_pressure
688 ! Define Clouds
689 if(imp_physics==99 .or. imp_physics==98)then
690 atmosphere(1)%cloud(1)%n_layers = lm
691 atmosphere(1)%cloud(1)%Type = water_cloud
692 atmosphere(1)%cloud(2)%n_layers = lm
693 atmosphere(1)%cloud(2)%Type = ice_cloud
694 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
695 atmosphere(1)%cloud(1)%n_layers = lm
696 atmosphere(1)%cloud(1)%Type = water_cloud
697 atmosphere(1)%cloud(2)%n_layers = lm
698 atmosphere(1)%cloud(2)%Type = ice_cloud
699 atmosphere(1)%cloud(3)%n_layers = lm
700 atmosphere(1)%cloud(3)%Type = rain_cloud
701 atmosphere(1)%cloud(4)%n_layers = lm
702 atmosphere(1)%cloud(4)%Type = snow_cloud
703 atmosphere(1)%cloud(5)%n_layers = lm
704 atmosphere(1)%cloud(5)%Type = graupel_cloud
705 atmosphere(1)%cloud(6)%n_layers = lm
706 atmosphere(1)%cloud(6)%Type = hail_cloud
707 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 .or. imp_physics==28 &
708 .or. imp_physics==11)then
709 atmosphere(1)%cloud(1)%n_layers = lm
710 atmosphere(1)%cloud(1)%Type = water_cloud
711 atmosphere(1)%cloud(2)%n_layers = lm
712 atmosphere(1)%cloud(2)%Type = ice_cloud
713 atmosphere(1)%cloud(3)%n_layers = lm
714 atmosphere(1)%cloud(3)%Type = rain_cloud
715 atmosphere(1)%cloud(4)%n_layers = lm
716 atmosphere(1)%cloud(4)%Type = snow_cloud
717 atmosphere(1)%cloud(5)%n_layers = lm
718 atmosphere(1)%cloud(5)%Type = graupel_cloud
719 end if
720
721 ! if(nchanl /= channelinfo(sensorindex)%n_channels) write(6,*)'***ERROR** nchanl,n_channels ', &
722 ! nchanl,channelinfo(sensorindex)%n_channels
723
724 ! Load surface sensor data structure
725 surface(1)%sensordata%n_channels = channelinfo(sensorindex)%n_channels
726 surface(1)%sensordata%sensor_id = channelinfo(sensorindex)%sensor_id
727 surface(1)%sensordata%WMO_sensor_id = channelinfo(sensorindex)%WMO_sensor_id
728 surface(1)%sensordata%WMO_Satellite_id = channelinfo(sensorindex)%WMO_Satellite_id
729 surface(1)%sensordata%sensor_channel = channelinfo(sensorindex)%sensor_channel
730
731 ! run crtm for nadir instruments / channels
732 nadir: if ( (isis=='imgr_g12' .and. (iget(327)>0 .or. &
733 iget(328)>0 .or. iget(329)>0 .or. iget(330)>0)) .or. &
734 (isis=='imgr_g11' .and. (iget(446)>0 .or. &
735 iget(447)>0 .or. iget(448)>0 .or. iget(449)>0)) .or. &
736 (isis=='amsre_aqua' .and. (iget(483) > 0 .or. iget(484) > 0 &
737 .or. iget(485) > 0 .or. iget(486) > 0)) .OR. &
738 (isis=='tmi_trmm' .and. (iget(488) > 0 .or. iget(489) > 0 &
739 .or. iget(490) > 0 .or. iget(491) > 0)) .OR. &
740 (isis=='abi_gr' .and. post_abigr) )then
741
742 do j=jsta,jend
743 loopi1:do i=ista,iend
744
745 ! Skiping the grids with filling value spval
746 do k=1,lm
747 if(abs(pmid(i,j,k)-spval)<=small .or. &
748 abs(t(i,j,k)-spval)<=small) then
749 do n=1,channelinfo(sensorindex)%n_channels
750 tb(i,j,n)=spval
751 enddo
752 cycle loopi1
753 endif
754 enddo
755
756 ! Load geometry structure
757 ! geometryinfo(1)%sensor_zenith_angle = zasat*rtd ! local zenith angle ???????
758 geometryinfo(1)%sensor_zenith_angle=0.
759 geometryinfo(1)%sensor_scan_angle=0.
760
761 !only call crtm if we have right satellite zenith angle
762 IF(geometryinfo(1)%sensor_zenith_angle <= max_sensor_scan_angle &
763 .and. geometryinfo(1)%sensor_zenith_angle >= 0.0_r_kind)THEN
764 geometryinfo(1)%source_zenith_angle = acos(czen(i,j))*rtd ! solar zenith angle
765 geometryinfo(1)%sensor_scan_angle = 0. ! scan angle, assuming nadir
766 if(i==ii.and.j==jj.and.debugprint)print*,'sample geometry ', &
767 geometryinfo(1)%sensor_zenith_angle &
768 ,geometryinfo(1)%source_zenith_angle &
769 ,czen(i,j)*rtd
770 ! Set land/sea, snow, ice percentages and flags
771 if(modelname == 'NCAR' .OR. modelname == 'RAPR')then
772 sfcpct(4)=pctsno(i,j)
773 else if(ivegsrc==1)then
774 itype=ivgtyp(i,j)
775 IF(itype == 0)itype=8
776 if(sno(i,j)<spval)then
777 snoeqv=sno(i,j)
778 else
779 snoeqv=0.
780 end if
781 CALL snfrac (snoeqv,ivgtyp(i,j),snofrac)
782 sfcpct(4)=snofrac
783 else if(ivegsrc==2)then
784 itype=ivgtyp(i,j)
785 itype = min(max(0,ivgtyp(i,j)),13)
786 if(sno(i,j)<spval)then
787 snoeqv=sno(i,j)
788 else
789 snoeqv=0.
790 end if
791 if(i==ii.and.j==jj.and.debugprint)print*,'sno,itype,ivgtyp B cing snfrc = ', &
792 snoeqv,itype,ivgtyp(i,j)
793 if(sm(i,j) > 0.1)then
794 sfcpct(4)=0.
795 else
796 call snfrac_gfs(snoeqv,ivgtyp(i,j),snofrac)
797 sfcpct(4)=snofrac
798 end if
799 end if
800
801! if (MODELNAME == 'GFS')then ! GFS uses 13 veg types
802! itype=IVGTYP(I,J)
803! itype = min(max(0,ivgtyp(i,j)),13)
804! ! IF(itype <= 0 .or. itype > 13)itype=7 !use scrub for ocean point
805! if(sno(i,j)/=spval)then
806! snoeqv=sno(i,j)
807! else
808! snoeqv=0.
809! end if
810! if(i==ii.and.j==jj)print*,'sno,itype,ivgtyp B cing snfrc = ', &
811! snoeqv,itype,IVGTYP(I,J)
812! if(sm(i,j) > 0.1)then
813! sfcpct(4)=0.
814! else
815! call snfrac_gfs(SNOeqv,IVGTYP(I,J),snofrac)
816! sfcpct(4)=snofrac
817! end if
818! if(i==ii.and.j==jj)print*,'sno,itype,ivgtyp,sfcpct(4) = ', &
819! snoeqv,itype,IVGTYP(I,J),sfcpct(4)
820! else if(MODELNAME == 'NCAR' .OR. MODELNAME == 'RAPR')then
821! sfcpct(4)=pctsno(i,j)
822! else
823! itype=IVGTYP(I,J)
824! IF(itype == 0)itype=8
825! CALL SNFRAC (SNO(I,J),IVGTYP(I,J),snofrac)
826! sfcpct(4)=snofrac
827! end if
828 ! CALL SNFRAC (SNO(I,J),IVGTYP(I,J),snofrac)
829 ! sfcpct(4)=snofrac
830 if(sm(i,j) > 0.1)then ! water
831 ! tsfc=sst(i,j)
832 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
833 vegcover=0.0
834 if(sfcpct(4) > 0.0_r_kind)then ! snow and water
835 sfcpct(1) = 1.0_r_kind-sfcpct(4)
836 sfcpct(2) = 0.0_r_kind
837 sfcpct(3) = 0.0_r_kind
838 else ! pure water
839 sfcpct(1) = 1.0_r_kind
840 sfcpct(2) = 0.0_r_kind
841 sfcpct(3) = 0.0_r_kind
842 end if
843 else ! land and sea ice
844 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
845 vegcover=vegfrc(i,j)
846 if(sice(i,j) > 0.1)then ! sea ice
847 if(sfcpct(4) > 0.0_r_kind)then ! sea ice and snow
848 sfcpct(3) = 1.0_r_kind-sfcpct(4)
849 sfcpct(1) = 0.0_r_kind
850 sfcpct(2) = 0.0_r_kind
851 else ! pure sea ice
852 sfcpct(3)= 1.0_r_kind
853 sfcpct(1) = 0.0_r_kind
854 sfcpct(2) = 0.0_r_kind
855 end if
856 else ! land
857 if(sfcpct(4) > 0.0_r_kind)then ! land and snow
858 sfcpct(2)= 1.0_r_kind-sfcpct(4)
859 sfcpct(1) = 0.0_r_kind
860 sfcpct(3) = 0.0_r_kind
861 else ! pure land
862 sfcpct(2)= 1.0_r_kind
863 sfcpct(1) = 0.0_r_kind
864 sfcpct(3) = 0.0_r_kind
865 end if
866 end if
867 end if
868 if(si(i,j)/=spval)then
869 snodepth = si(i,j)
870 else
871 snodepth = 0.
872 end if
873 ! Chuang: for igbp type 15 (snow/ice), the main type needs to be set to ice or snow
874 ! to prevent crtm forward model from failing
875 if(ivegsrc==1 .and. itype==15 .and. sfcpct(4)<1.0_r_kind)then
876 if(debugprint)print*,'changing land type for veg type 15',i,j,itype,sfcpct(1:4)
877 sfcpct(1)=0.0_r_kind
878 sfcpct(2)=0.0_r_kind
879 sfcpct(3)=0.0_r_kind
880 sfcpct(4)=1.0_r_kind
881 !print*,'change main land type to snow for veg type 15 ',i,j
882 end if
883
884 sea = sfcpct(1) >= 0.99_r_kind
885 land = sfcpct(2) >= 0.99_r_kind
886 ice = sfcpct(3) >= 0.99_r_kind
887 snow = sfcpct(4) >= 0.99_r_kind
888 mixed = .not. sea .and. .not. ice .and. &
889 .not. land .and. .not. snow
890 if((sfcpct(1)+sfcpct(2)+sfcpct(3)+sfcpct(4)) >1._r_kind) &
891 print*,'ERROR sfcpct ',i,j,sfcpct(1) &
892 ,sfcpct(2),sfcpct(3),sfcpct(4)
893 ! Load surface structure
894
895 ! Define land characteristics
896
897 ! **NOTE: The model surface type --> CRTM surface type
898 ! mapping below is specific to the versions NCEP
899 ! GFS and NNM as of September 2005
900 ! itype = ivgtyp(i,j)
901 if(ivegsrc==0)then
902 itype = min(max(0,ivgtyp(i,j)),novegtype)
903 else
904 itype = min(max(1,ivgtyp(i,j)),novegtype)
905 end if
906 surface(1)%land_type = model_to_crtm(itype)
907
908 if(gridtype=='B' .or. gridtype=='E')then
909 surface(1)%wind_speed = sqrt(u10h(i,j)*u10h(i,j) &
910 +v10h(i,j)*v10h(i,j))
911 else
912 surface(1)%wind_speed = sqrt(u10(i,j)*u10(i,j) &
913 +v10(i,j)*v10(i,j))
914 end if
915 surface(1)%water_coverage = sfcpct(1)
916 surface(1)%land_coverage = sfcpct(2)
917 surface(1)%ice_coverage = sfcpct(3)
918 surface(1)%snow_coverage = sfcpct(4)
919
920 surface(1)%land_temperature = tsfc
921 surface(1)%snow_temperature = min(tsfc,280._r_kind)
922 surface(1)%water_temperature = max(tsfc,270._r_kind)
923 surface(1)%ice_temperature = min(tsfc,280._r_kind)
924 if(smstot(i,j)/=spval)then
925 surface(1)%soil_moisture_content = smstot(i,j)/10. !convert to cgs !???
926 else
927 surface(1)%soil_moisture_content = 0.05 ! default crtm value
928 end if
929 surface(1)%vegetation_fraction = vegcover
930 ! surface(1)%vegetation_fraction = vegfrc(i,j)
931 surface(1)%soil_temperature = 283.
932 ! surface(1)%soil_temperature = stc(i,j,1)
933 surface(1)%snow_depth = snodepth ! in mm
934 ! Debug print
935 if(debugprint)then
936 if(surface(1)%wind_speed<0. .or. surface(1)%wind_speed>200.) &
937 print*,'bad 10 m wind'
938 if(surface(1)%water_coverage<0. .or. surface(1)%water_coverage>1.) &
939 print*,'bad water coverage'
940 if(surface(1)%land_coverage<0. .or. surface(1)%land_coverage>1.) &
941 print*,'bad land coverage'
942 if(surface(1)%ice_coverage<0. .or. surface(1)%ice_coverage>1.) &
943 print*,'bad ice coverage'
944 if(surface(1)%snow_coverage<0. .or. surface(1)%snow_coverage>1.) &
945 print*,'bad snow coverage'
946 if(surface(1)%land_temperature<0. .or. surface(1)%land_temperature>350.) &
947 print*,'bad land T'
948 if(surface(1)%soil_moisture_content<0. .or. surface(1)%soil_moisture_content>600.) &
949 print*,'bad soil_moisture_content'
950 if(surface(1)%vegetation_fraction<0. .or. surface(1)%vegetation_fraction>1.) &
951 print*,'bad vegetation cover'
952 if(surface(1)%snow_depth<0. .or. surface(1)%snow_depth>10000.) &
953 print*,'bad snow_depth'
954 end if
955 if(i==ii.and.j==jj.and.debugprint)print*,'sample surface in CALRAD=', &
956 i,j,surface(1)%wind_speed,surface(1)%water_coverage, &
957 surface(1)%land_coverage,surface(1)%ice_coverage, &
958 surface(1)%snow_coverage,surface(1)%land_temperature, &
959 surface(1)%snow_temperature,surface(1)%water_temperature, &
960 surface(1)%ice_temperature,surface(1)%vegetation_fraction, &
961 surface(1)%soil_temperature,surface(1)%snow_depth, &
962 surface(1)%land_type,sm(i,j)
963
964 ! Load profiles into model layers
965
966 ! Load atmosphere profiles into RTM model layers
967 ! CRTM counts from top down just as post does
968 if(i==ii.and.j==jj.and.debugprint)print*,'TOA= ',atmosphere(1)%level_pressure(0)
969 do k = 1,lm
970 atmosphere(1)%cloud_fraction(k) = min(max(cfr(i,j,k),0.),1.)
971 atmosphere(1)%level_pressure(k) = pint(i,j,k+1)/r100
972 atmosphere(1)%pressure(k) = pmid(i,j,k)/r100
973 atmosphere(1)%temperature(k) = t(i,j,k)
974 atmosphere(1)%absorber(k,1) = max(0. &
975 ,q(i,j,k)*h1000/(h1-q(i,j,k))) ! use mixing ratio like GSI
976 atmosphere(1)%absorber(k,2) = max(ozsmall,o3(i,j,k)*constoz)
977 ! atmosphere(1)%absorber(k,2) = max(ozsmall,o3(i,j,k)*h1000) ! convert to g/kg
978 ! fill in cloud mixing ratio later
979 if(debugprint)then
980 if(atmosphere(1)%level_pressure(k)<0. .or. atmosphere(1)%level_pressure(k)>1060.) &
981 & print*,'bad atmosphere(1)%level_pressure' &
982 & ,i,j,k,atmosphere(1)%level_pressure(k)
983 if(atmosphere(1)%pressure(k)<0. .or. &
984 & atmosphere(1)%pressure(k)>1060.) &
985 & print*,'bad atmosphere(1)%pressure' &
986 & ,i,j,k,atmosphere(1)%pressure(k)
987 if(atmosphere(1)%temperature(k)<0. .or. &
988 & atmosphere(1)%temperature(k)>400.) &
989 & print*,'bad atmosphere(1)%temperature'
990 ! if(atmosphere(1)%absorber(k,1)<0. .or. &
991 ! & atmosphere(1)%absorber(k,1)>1.) &
992 ! & print*,'bad atmosphere water vapor'
993 ! if(atmosphere(1)%absorber(k,2)<0. .or. &
994 ! & atmosphere(1)%absorber(k,1)>1.) &
995 ! & print*,'bad atmosphere o3'
996 end if
997! if(i==ii.and.j==jj.and.debugprint)print*,'sample atmosphere in CALRAD=', &
998! i,j,k,atmosphere(1)%level_pressure(k),atmosphere(1)%pressure(k), &
999! atmosphere(1)%temperature(k),atmosphere(1)%absorber(k,1), &
1000! atmosphere(1)%absorber(k,2)
1001 ! Specify clouds
1002 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g !crtm uses column integrated field
1003 if(imp_physics==99 .or. imp_physics==98)then
1004 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1005 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1006 ! GFS uses temperature and ice concentration dependency formulation to
1007 ! determine effetive radis for cloud ice
1008 ! since GFS does not output ice concentration yet, use default 50 um
1009 atmosphere(1)%cloud(2)%effective_radius(k) = 50.
1010 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1011 if(debugprint)then
1012 if(atmosphere(1)%cloud(1)%water_content(k)<0. .or. &
1013 atmosphere(1)%cloud(1)%water_content(k)>1.) &
1014 print*,'bad atmosphere cloud water'
1015 if(atmosphere(1)%cloud(2)%water_content(k)<0. .or. &
1016 atmosphere(1)%cloud(2)%water_content(k)>1.) &
1017 print*,'bad atmosphere cloud ice'
1018 end if
1019 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
1020 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1021 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1022 atmosphere(1)%cloud(2)%effective_radius(k) = 75.
1023 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1024 rhox=1000.
1025 rho=pmid(i,j,k)/(rd*t(i,j,k)*(1.+d608*q(i,j,k)))
1026 atmosphere(1)%cloud(3)%effective_radius(k) = 0.
1027 if(nrain(i,j,k)>0.) &
1028 atmosphere(1)%cloud(3)%effective_radius(k) = &
1029 1.0e6*1.5*(rho*qqr(i,j,k)/(pi*rhox*nrain(i,j,k)))**(1./3.)
1030 atmosphere(1)%cloud(3)%water_content(k) = max(0.,qqr(i,j,k)*dpovg)
1031 if(f_rimef(i,j,k)<=5.0)then
1032 rhox=100
1033 if(nlice(i,j,k)>0.) &
1034 atmosphere(1)%cloud(4)%effective_radius(k) = &
1035 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.) !convert to microns
1036 atmosphere(1)%cloud(4)%water_content(k) = max(0.,qqs(i,j,k)*dpovg)
1037 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1038 atmosphere(1)%cloud(5)%water_content(k) =0.
1039 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1040 atmosphere(1)%cloud(6)%water_content(k) =0.
1041 else if(f_rimef(i,j,k)<=20.0)then
1042 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1043 atmosphere(1)%cloud(4)%water_content(k) =0.
1044 rhox=400.
1045 if(nlice(i,j,k)>0.) &
1046 atmosphere(1)%cloud(5)%effective_radius(k) = &
1047 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1048 atmosphere(1)%cloud(5)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1049 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1050 atmosphere(1)%cloud(6)%water_content(k) =0.
1051 else
1052 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1053 atmosphere(1)%cloud(4)%water_content(k) =0.
1054 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1055 atmosphere(1)%cloud(5)%water_content(k) =0.
1056 rhox=900.
1057 if(nlice(i,j,k)>0.) &
1058 atmosphere(1)%cloud(6)%effective_radius(k) = &
1059 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1060 atmosphere(1)%cloud(6)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1061 end if
1062 if(debugprint .and. i==im/2 .and. j==jsta) &
1063 print*,'sample precip ice radius= ',i,j,k, f_rimef(i,j,k), &
1064 atmosphere(1)%cloud(4)%effective_radius(k), atmosphere(1)%cloud(4)%water_content(k), &
1065 atmosphere(1)%cloud(5)%effective_radius(k), atmosphere(1)%cloud(5)%water_content(k), &
1066 atmosphere(1)%cloud(6)%effective_radius(k), atmosphere(1)%cloud(6)%water_content(k)
1067
1068 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 .or. &
1069 imp_physics==28 .or. imp_physics==11)then
1070 atmosphere(1)%cloud(1)%water_content(k)=max(0.,qqw(i,j,k)*dpovg)
1071 atmosphere(1)%cloud(2)%water_content(k)=max(0.,qqi(i,j,k)*dpovg)
1072 atmosphere(1)%cloud(3)%water_content(k)=max(0.,qqr(i,j,k)*dpovg)
1073 atmosphere(1)%cloud(4)%water_content(k)=max(0.,qqs(i,j,k)*dpovg)
1074 atmosphere(1)%cloud(5)%water_content(k)=max(0.,qqg(i,j,k)*dpovg)
1075 atmosphere(1)%cloud(1)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1076 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1077 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'C')
1078 atmosphere(1)%cloud(2)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1079 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1080 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'I')
1081 atmosphere(1)%cloud(3)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1082 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1083 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'R')
1084 atmosphere(1)%cloud(4)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1085 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1086 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'S')
1087 atmosphere(1)%cloud(5)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1088 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1089 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'G')
1090 end if
1091 end do
1092!Meng 09/2018 modify two model layer having identical pressure
1093 do k = 1,lm-1
1094 if (abs(atmosphere(1)%pressure(k)-atmosphere(1)%pressure(k+1)) &
1095 < 0.005) then
1096 atmosphere(1)%pressure(k)=atmosphere(1)%pressure(k)+0.005
1097 endif
1098 enddo
1099
1100 !bsf - start
1101 !-- Add subgrid-scale convective clouds for WRF runs
1102 if(icu_physics==2) then
1103 lcbot=nint(hbot(i,j))
1104 lctop=nint(htop(i,j))
1105 if (lcbot-lctop > 1) then
1106 !-- q_conv = assumed grid-averaged cloud water/ice condensate from conimp_physics,vection (Cu)
1107 ! In "params" Qconv=0.1e-3 and TRAD_ice=253.15; cnvcfr is the Cu cloud fraction
1108 ! calculated as a function of Cu rain rate (Slingo, 1987) in subroutine MDLFLD
1109 q_conv = cnvcfr(i,j)*qconv
1110 do k = lctop,lcbot
1111 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g
1112 if (t(i,j,k) < trad_ice) then
1113 atmosphere(1)%cloud(2)%water_content(k) = &
1114 atmosphere(1)%cloud(2)%water_content(k) + dpovg*q_conv
1115 else
1116 atmosphere(1)%cloud(1)%water_content(k) = &
1117 atmosphere(1)%cloud(1)%water_content(k) + dpovg*q_conv
1118 endif
1119 end do !-- do k = lctop,lcbot
1120 endif !-- if (lcbot-lctop > 1) then
1121 endif !-- if (MODELNAME == 'NMM' .OR. MODELNAME == 'NCAR') then
1122 !bsf - end
1123
1124 ! call crtm forward model
1125 error_status = crtm_forward(atmosphere,surface, &
1126 geometryinfo,channelinfo(sensorindex:sensorindex), &
1127 rtsolution)
1128 if (error_status /=0) then
1129 print*,'***ERROR*** during crtm_forward call ', error_status
1130 do n=1,channelinfo(sensorindex)%n_channels
1131 tb(i,j,n)=spval
1132 end do
1133 ! tb2(i,j)=spval
1134 ! tb3(i,j)=spval
1135 ! tb4(i,j)=spval
1136 else
1137 do n=1,channelinfo(sensorindex)%n_channels
1138 tb(i,j,n)=rtsolution(n,1)%brightness_temperature
1139 end do
1140 ! tb1(i,j)=rtsolution(1,1)%brightness_temperature
1141 ! tb2(i,j)=rtsolution(2,1)%brightness_temperature
1142 ! tb3(i,j)=rtsolution(3,1)%brightness_temperature
1143 ! tb4(i,j)=rtsolution(4,1)%brightness_temperature
1144 if(i==ii.and.j==jj) then
1145 do n=1,channelinfo(sensorindex)%n_channels
1146 3301 format('Sample rtsolution(',i0,',',i0,') in CALRAD = ',f0.3)
1147 print 3301,n,1,rtsolution(n,1)%brightness_temperature
1148 enddo
1149 do n=1,channelinfo(sensorindex)%n_channels
1150 3302 format('Sample tb(',i0,',',i0,',',i0,') in CALRAD = ',f0.3)
1151 print 3302,ii,jj,n,tb(ii,jj,n)
1152 enddo
1153 endif
1154 ! if(tb1(i,j) < 400. ) &
1155 ! & print*,'good tb1 ',i,j,tb1(i,j),gdlat(i,j),gdlon(i,j)
1156 ! if(tb2(i,j) > 400.)print*,'bad tb2 ',i,j,tb2(i,j)
1157 ! if(tb3(i,j) > 400.)print*,'bad tb3 ',i,j,tb3(i,j)
1158 ! if(tb4(i,j) > 400.)print*,'bad tb4 ',i,j,tb4(i,j)
1159 end if
1160 else
1161 do n=1,channelinfo(sensorindex)%n_channels
1162 tb(i,j,n)=spval
1163 end do
1164 ! tb1(i,j)=spval
1165 ! tb2(i,j)=spval
1166 ! tb3(i,j)=spval
1167 ! tb4(i,j)=spval
1168 END IF ! endif block for allowable satellite zenith angle
1169 end do loopi1 ! end loop for i
1170 end do ! end loop for j
1171
1172 ! error_status = crtm_destroy(channelinfo)
1173 ! if (error_status /= success) &
1174 ! & print*,'ERROR*** crtm_destroy error_status=',error_status
1175
1176 if (isis=='amsre_aqua')then ! writing amsre to grib (37 & 89 GHz)
1177 do ixchan=1,4
1178 ichan=8+ixchan
1179 igot=iget(482+ixchan)
1180 if(igot>0) then
1181 do j=jsta,jend
1182 do i=ista,iend
1183 grid1(i,j)=tb(i,j,ichan)
1184 enddo
1185 enddo
1186 if (grib=="grib2") then
1187 cfld=cfld+1
1188 fld_info(cfld)%ifld=iavblfld(igot)
1189 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1190 endif
1191 endif
1192 enddo
1193 end if ! end of outputting amsre
1194 if (isis=='tmi_trmm')then ! writing trmm to grib (37 & 85.5 GHz)
1195 do ixchan=1,4
1196 ichan=5+ixchan
1197 igot=iget(487+ixchan)
1198 if(igot>0) then
1199 do j=jsta,jend
1200 do i=ista,iend
1201 grid1(i,j) = tb(i,j,ichan)
1202 enddo
1203 enddo
1204 if (grib=="grib2") then
1205 cfld=cfld+1
1206 fld_info(cfld)%ifld=iavblfld(igot)
1207 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1208 endif
1209 endif
1210 enddo
1211 end if ! end of outputting trmm
1212
1213 if (isis=='imgr_g11')then ! writing goes 11 to grib
1214 do ixchan=1,4
1215 ichan=ixchan
1216 igot=445+ixchan
1217 if(igot>0) then
1218 do j=jsta,jend
1219 do i=ista,iend
1220 grid1(i,j) = tb(i,j,ichan)
1221 enddo
1222 enddo
1223 if (grib=="grib2") then
1224 cfld=cfld+1
1225 fld_info(cfld)%ifld=iavblfld(igot)
1226 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1227 endif
1228 endif ! IGOT
1229 enddo
1230 end if ! end of outputting goes 11
1231
1232 if (isis=='imgr_g12')then ! writing goes 12 to grib
1233 do ixchan=1,4 ! write brightness temperatures
1234 ichan=ixchan
1235 igot=iget(326+ixchan)
1236 if(igot>0) then
1237 do j=jsta,jend
1238 do i=ista,iend
1239 grid1(i,j)=tb(i,j,ichan)
1240 enddo
1241 enddo
1242 if (grib=="grib2") then
1243 cfld=cfld+1
1244 fld_info(cfld)%ifld=iavblfld(igot)
1245 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1246 endif
1247 endif
1248 enddo
1249 endif ! end of outputting goes 12
1250 if (isis=='abi_gr')then ! writing goes-r nadir to grib2
1251 nc=0
1252 do ixchan=1,10
1253 ichan=ixchan
1254 igot=iget(957+ixchan)
1255 if(igot>0)then
1256 do j=jsta,jend
1257 do i=ista,iend
1258 grid1(i,j)=tb(i,j,ichan)
1259 enddo
1260 enddo
1261 if(grib=="grib2" )then
1262 cfld=cfld+1
1263 fld_info(cfld)%ifld=iavblfld(igot)
1264 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1265 endif
1266 endif
1267 enddo ! channel loop
1268 end if ! end of outputting goes-r nadir
1269
1270
1271 end if nadir ! end if for computing nadir simulated radiance
1272
1273 ! run crtm for non-nadir instruments / channels
1274
1275 nonnadir: if((isis=='ssmi_f13' .and. iget(800) > 0 ) .OR. &
1276 (isis=='ssmi_f14' .and. iget(806) > 0 ) .OR. &
1277 (isis=='ssmi_f15' .and. iget(812) > 0 ) .OR. &
1278 (isis=='ssmis_f16' .and. iget(818) > 0) .OR. &
1279 (isis=='ssmis_f17' .and. post_ssmis17) .OR. &
1280 (isis=='ssmis_f18' .and. iget(832) > 0) .OR. &
1281 (isis=='ssmis_f19' .and. iget(839) > 0) .OR. &
1282 (isis=='ssmis_f20' .and. iget(846) > 0) .OR. &
1283 (isis=='imgr_mt2' .and. iget(860)>0) .OR. &
1284 (isis=='imgr_mt1r' .and. iget(864)>0) .OR. &
1285 (isis=='imgr_insat3d' .and. iget(865)>0) .OR. &
1286 (isis=='imgr_g13' .and. iget(868)>0) .OR. &
1287 (isis=='imgr_g15' .and. iget(872)>0) .OR. &
1288 (isis=='abi_g16' .and. post_abig16) .OR. &
1289 (isis=='abi_g17' .and. post_abig17) .OR. &
1290 (isis=='abi_g18' .and. post_abig18) .OR. &
1291 (isis=='seviri_m10' .and. iget(876)>0) .OR. &
1292 (isis=='ahi_himawari8' .and. post_ahi8) .OR. &
1293 (isis=='imgr_g12' .and. (iget(456)>0 .or. &
1294 iget(457)>0 .or. iget(458)>0 .or. iget(459)>0)) .or. &
1295 (isis=='imgr_g11' .and. (iget(460)>0 .or. &
1296 iget(461)>0 .or. iget(462)>0 .or. iget(463)>0)))then
1297
1298 do j=jsta,jend
1299 loopi2:do i=ista,iend
1300
1301 ! Skiping the grids with filling value spval
1302 do k=1,lm
1303 if(abs(pmid(i,j,k)-spval)<=small .or. &
1304 abs(t(i,j,k)-spval)<=small) then
1305 do n=1,channelinfo(sensorindex)%n_channels
1306 tb(i,j,n)=spval
1307 enddo
1308 cycle loopi2
1309 endif
1310 enddo
1311
1312 ! Load geometry structure
1313 ! geometryinfo(1)%sensor_zenith_angle = zasat*rtd ! local zenith angle ???????
1314 ! compute satellite zenith angle
1315 if(isis=='imgr_g12')then
1316 sublat=0.0
1317 sublon=-75.0
1318 else if(isis=='seviri_m10')then
1319 sublat=0.0
1320 sublon=0.0
1321 else if(isis=='imgr_g13')then
1322 sublat=0.0
1323 sublon=-75.0
1324 else if(isis=='imgr_g15')then
1325 sublat=0.0
1326 sublon=-135.0
1327 else if(isis=='abi_g16')then ! positions should be controlled by runtime setting or fix file
1328 sublat=0.0
1329 sublon=-75.2
1330 else if(isis=='abi_g17')then
1331 sublat=0.0
1332 sublon=-137.2
1333 else if(isis=='abi_g18')then
1334 sublat=0.0
1335 sublon=-137.2
1336 else if(isis=='imgr_g11')then
1337 sublat=0.0
1338 sublon=-135.0
1339 else if(isis=='imgr_mt2') then
1340 sublat=0.0
1341 sublon=145.0
1342 else if(isis=='imgr_mt1r') then
1343 sublat=0.0
1344 sublon=140.0
1345 else if(isis=='imgr_insat3d') then
1346 sublat=0.0
1347 sublon=74.0
1348 else if(isis=='ahi_himawari8') then
1349 sublat=0.0
1350 sublon=140.7
1351 end if
1352
1353! use zenith angle = 53.1 for SSMI and SSMIS:
1354 if(isis=='ssmis_f16'.or.isis=='ssmis_f17'.or.isis=='ssmis_f18'.or. &
1355 isis=='ssmis_f19'.or.isis=='ssmis_f20'.or.isis=='ssmi_f13'.or. &
1356 isis=='ssmi_f14'.or.isis=='ssmi_f15')then
1357 sat_zenith=53.1
1358 else
1359 ! For other imagers (GOES-11 and 12), calculate based on satellite location:
1360 call geo_zenith_angle(i,j,gdlat(i,j),gdlon(i,j) &
1361 ,sublat,sublon,sat_zenith)
1362 endif
1363
1364 geometryinfo(1)%sensor_zenith_angle=sat_zenith
1365 geometryinfo(1)%sensor_scan_angle=sat_zenith
1366
1367 if(i==ii .and. j==jj.and.debugprint) then
1368 print *,'zenith info: zenith=',sat_zenith,' scan=',sat_zenith, &
1369 ' MAX_SENSOR_SCAN_ANGLE=',max_sensor_scan_angle
1370 endif
1371 ! geometryinfo(1)%sensor_zenith_angle = 0. ! 44.
1372 !only call crtm if we have right satellite zenith angle
1373 IF(geometryinfo(1)%sensor_zenith_angle <= max_sensor_scan_angle &
1374 .and. geometryinfo(1)%sensor_zenith_angle >= 0.0_r_kind)THEN
1375 geometryinfo(1)%source_zenith_angle = acos(czen(i,j))*rtd ! solar zenith angle
1376 geometryinfo(1)%sensor_scan_angle = 0. ! scan angle, assuming nadir
1377 if(i==ii.and.j==jj.and.debugprint)print*,'sample geometry ', &
1378 geometryinfo(1)%sensor_zenith_angle &
1379 ,geometryinfo(1)%source_zenith_angle &
1380 ,czen(i,j)*rtd
1381 ! Set land/sea, snow, ice percentages and flags
1382
1383 if(modelname == 'NCAR' .OR. modelname == 'RAPR')then
1384 sfcpct(4)=pctsno(i,j)
1385 else if(ivegsrc==1)then
1386 itype=ivgtyp(i,j)
1387 IF(itype == 0)itype=8
1388 if(sno(i,j)<spval)then
1389 snoeqv=sno(i,j)
1390 else
1391 snoeqv=0.
1392 end if
1393 CALL snfrac (snoeqv,ivgtyp(i,j),snofrac)
1394 sfcpct(4)=snofrac
1395 else if(ivegsrc==2)then
1396 itype=ivgtyp(i,j)
1397 itype = min(max(0,ivgtyp(i,j)),13)
1398 if(sno(i,j)<spval)then
1399 snoeqv=sno(i,j)
1400 else
1401 snoeqv=0.
1402 end if
1403 if(i==ii.and.j==jj)print*,'sno,itype,ivgtyp B cing snfrc = ', &
1404 snoeqv,itype,ivgtyp(i,j)
1405 if(sm(i,j) > 0.1)then
1406 sfcpct(4)=0.
1407 else
1408 call snfrac_gfs(snoeqv,ivgtyp(i,j),snofrac)
1409 sfcpct(4)=snofrac
1410 end if
1411 end if
1412
1413 ! CALL SNFRAC (SNO(I,J),IVGTYP(I,J),snofrac)
1414 ! sfcpct(4)=snofrac
1415 if(sm(i,j) > 0.1)then ! water
1416 ! tsfc=sst(i,j)
1417 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
1418 vegcover=0.0
1419 if(sfcpct(4) > 0.0_r_kind)then ! snow and water
1420 sfcpct(1) = 1.0_r_kind-sfcpct(4)
1421 sfcpct(2) = 0.0_r_kind
1422 sfcpct(3) = 0.0_r_kind
1423 else ! pure water
1424 sfcpct(1) = 1.0_r_kind
1425 sfcpct(2) = 0.0_r_kind
1426 sfcpct(3) = 0.0_r_kind
1427 end if
1428 else ! land and sea ice
1429 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
1430 vegcover=vegfrc(i,j)
1431 if(sice(i,j) > 0.1)then ! sea ice
1432 if(sfcpct(4) > 0.0_r_kind)then ! sea ice and snow
1433 sfcpct(3) = 1.0_r_kind-sfcpct(4)
1434 sfcpct(1) = 0.0_r_kind
1435 sfcpct(2) = 0.0_r_kind
1436 else ! pure sea ice
1437 sfcpct(3)= 1.0_r_kind
1438 sfcpct(1) = 0.0_r_kind
1439 sfcpct(2) = 0.0_r_kind
1440 end if
1441 else ! land
1442 if(sfcpct(4) > 0.0_r_kind)then ! land and snow
1443 sfcpct(2)= 1.0_r_kind-sfcpct(4)
1444 sfcpct(1) = 0.0_r_kind
1445 sfcpct(3) = 0.0_r_kind
1446 else ! pure land
1447 sfcpct(2)= 1.0_r_kind
1448 sfcpct(1) = 0.0_r_kind
1449 sfcpct(3) = 0.0_r_kind
1450 end if
1451 end if
1452 end if
1453 if(si(i,j)/=spval)then
1454 snodepth = si(i,j)
1455 else
1456 snodepth = 0.
1457 end if
1458 !DTC added based on nadir section
1459 ! Chuang: for igbp type 15 (snow/ice), the main type
1460 ! needs to be set to ice or snow
1461 ! to prevent crtm forward model from failing
1462 if(novegtype==20 .and. itype==15 .and.sfcpct(4)<1.0_r_kind)then
1463 if(debugprint)print*,'changing land type for veg type 15',i,j,itype,sfcpct(1:4)
1464 sfcpct(1)=0.0_r_kind
1465 sfcpct(2)=0.0_r_kind
1466 sfcpct(3)=0.0_r_kind
1467 sfcpct(4)=1.0_r_kind
1468 !print*,'change main land type to snow for veg type 15
1469 !',i,j
1470 end if
1471
1472 sea = sfcpct(1) >= 0.99_r_kind
1473 land = sfcpct(2) >= 0.99_r_kind
1474 ice = sfcpct(3) >= 0.99_r_kind
1475 snow = sfcpct(4) >= 0.99_r_kind
1476 mixed = .not. sea .and. .not. ice .and. &
1477 .not. land .and. .not. snow
1478 if((sfcpct(1)+sfcpct(2)+sfcpct(3)+sfcpct(4)) >1._r_kind) &
1479 print*,'ERROR sfcpct ',i,j,sfcpct(1) &
1480 ,sfcpct(2),sfcpct(3),sfcpct(4)
1481 ! Load surface structure
1482
1483 ! Define land characteristics
1484
1485 ! **NOTE: The model surface type --> CRTM surface type
1486 ! mapping below is specific to the versions NCEP
1487 ! GFS and NNM as of September 2005
1488 ! itype = ivgtyp(i,j)
1489 if(ivegsrc==0)then
1490 itype = min(max(0,ivgtyp(i,j)),novegtype)
1491 else
1492 itype = min(max(1,ivgtyp(i,j)),novegtype)
1493 end if
1494 surface(1)%land_type = model_to_crtm(itype)
1495
1496 if(gridtype=='B' .or. gridtype=='E')then
1497 surface(1)%wind_speed = sqrt(u10h(i,j)*u10h(i,j) &
1498 +v10h(i,j)*v10h(i,j))
1499 else
1500 surface(1)%wind_speed = sqrt(u10(i,j)*u10(i,j) &
1501 +v10(i,j)*v10(i,j))
1502 end if
1503
1504 surface(1)%water_coverage = sfcpct(1)
1505 surface(1)%land_coverage = sfcpct(2)
1506 surface(1)%ice_coverage = sfcpct(3)
1507 surface(1)%snow_coverage = sfcpct(4)
1508 surface(1)%land_temperature = tsfc
1509 surface(1)%snow_temperature = min(tsfc,280._r_kind)
1510 surface(1)%water_temperature = max(tsfc,270._r_kind)
1511 surface(1)%ice_temperature = min(tsfc,280._r_kind)
1512 if(smstot(i,j)/=spval)then
1513 surface(1)%soil_moisture_content = smstot(i,j)/10. !convert to cgs !???
1514 else
1515 surface(1)%soil_moisture_content = 0.05 ! default crtm value
1516 end if
1517 surface(1)%vegetation_fraction = vegcover
1518 ! surface(1)%vegetation_fraction = vegfrc(i,j)
1519 surface(1)%soil_temperature = 283.
1520 ! surface(1)%soil_temperature = stc(i,j,1)
1521 surface(1)%snow_depth = snodepth ! in mm
1522 ! Debug print
1523 if(debugprint)then
1524 if(surface(1)%wind_speed<0. .or. surface(1)%wind_speed>200.) &
1525 print*,'bad 10 m wind'
1526 if(surface(1)%water_coverage<0. .or. surface(1)%water_coverage>1.) &
1527 print*,'bad water coverage'
1528 if(surface(1)%land_coverage<0. .or. surface(1)%land_coverage>1.) &
1529 print*,'bad land coverage'
1530 if(surface(1)%ice_coverage<0. .or. surface(1)%ice_coverage>1.) &
1531 print*,'bad ice coverage'
1532 if(surface(1)%snow_coverage<0. .or. surface(1)%snow_coverage>1.) &
1533 print*,'bad snow coverage'
1534 if(surface(1)%land_temperature<0. .or. surface(1)%land_temperature>350.) &
1535 print*,'bad land T'
1536 if(surface(1)%soil_moisture_content<0. .or. surface(1)%soil_moisture_content>600.) &
1537 print*,'bad soil_moisture_content'
1538 if(surface(1)%vegetation_fraction<0. .or. surface(1)%vegetation_fraction>1.) &
1539 print*,'bad vegetation cover'
1540 if(surface(1)%snow_depth<0. .or. surface(1)%snow_depth>10000.) &
1541 print*,'bad snow_depth'
1542 end if
1543
1544 if(i==ii.and.j==jj.and.debugprint)print*,'sample surface in CALRAD=', &
1545 i,j,surface(1)%wind_speed,surface(1)%water_coverage, &
1546 surface(1)%land_coverage,surface(1)%ice_coverage, &
1547 surface(1)%snow_coverage,surface(1)%land_temperature, &
1548 surface(1)%snow_temperature,surface(1)%water_temperature, &
1549 surface(1)%ice_temperature,surface(1)%vegetation_fraction, &
1550 surface(1)%soil_temperature,surface(1)%snow_depth, &
1551 surface(1)%land_type,sm(i,j)
1552
1553 ! Load profiles into model layers
1554
1555 ! Load atmosphere profiles into RTM model layers
1556 ! CRTM counts from top down just as post does
1557 if(i==ii.and.j==jj.and.debugprint)print*,'TOA= ',atmosphere(1)%level_pressure(0)
1558 do k = 1,lm
1559 atmosphere(1)%cloud_fraction(k) = min(max(cfr(i,j,k),0.),1.)
1560 atmosphere(1)%level_pressure(k) = pint(i,j,k+1)/r100
1561 atmosphere(1)%pressure(k) = pmid(i,j,k)/r100
1562 atmosphere(1)%temperature(k) = t(i,j,k)
1563 atmosphere(1)%absorber(k,1) = max(0. &
1564 ,q(i,j,k)*h1000/(h1-q(i,j,k))) ! use mixing ratio like GSI
1565 atmosphere(1)%absorber(k,2) = max(ozsmall,o3(i,j,k)*constoz)
1566 ! atmosphere(1)%absorber(k,2) = max(ozsmall,o3(i,j,k)*h1000) ! convert to g/kg
1567 ! fill in cloud mixing ratio later
1568 if(debugprint)then
1569 if(atmosphere(1)%level_pressure(k)<0. .or. atmosphere(1)%level_pressure(k)>1060.) &
1570 print*,'bad atmosphere(1)%level_pressure' &
1571 ,i,j,k,atmosphere(1)%level_pressure(k)
1572 if(atmosphere(1)%pressure(k)<0. .or. &
1573 atmosphere(1)%pressure(k)>1060.) &
1574 print*,'bad atmosphere(1)%pressure' &
1575 ,i,j,k,atmosphere(1)%pressure(k)
1576 if(atmosphere(1)%temperature(k)<0. .or. &
1577 atmosphere(1)%temperature(k)>400.) &
1578 print*,'bad atmosphere(1)%temperature'
1579 ! if(atmosphere(1)%absorber(k,1)<0. .or. &
1580 ! & atmosphere(1)%absorber(k,1)>1.) &
1581 ! & print*,'bad atmosphere water vapor'
1582 ! if(atmosphere(1)%absorber(k,2)<0. .or. &
1583 ! & atmosphere(1)%absorber(k,1)>1.) &
1584 ! & print*,'bad atmosphere o3'
1585 end if
1586! if(i==ii.and.j==jj)print*,'sample atmosphere in CALRAD=', &
1587! i,j,k,atmosphere(1)%level_pressure(k),atmosphere(1)%pressure(k), &
1588! atmosphere(1)%temperature(k),atmosphere(1)%absorber(k,1), &
1589! atmosphere(1)%absorber(k,2)
1590 ! Specify clouds
1591 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g !crtm uses column integrated field
1592 if(imp_physics==99 .or. imp_physics==98)then
1593 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1594 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1595 ! GFS uses temperature and ice concentration dependency formulation to determine
1596 ! effetive radis for cloud ice since GFS does not output ice concentration yet,
1597 !use default 50 um
1598 atmosphere(1)%cloud(2)%effective_radius(k) = 50.
1599 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1600 if(debugprint)then
1601 if(atmosphere(1)%cloud(1)%water_content(k)<0. .or. &
1602 atmosphere(1)%cloud(1)%water_content(k)>1.) &
1603 print*,'bad atmosphere cloud water'
1604 if(atmosphere(1)%cloud(2)%water_content(k)<0. .or. &
1605 atmosphere(1)%cloud(2)%water_content(k)>1.) &
1606 print*,'bad atmosphere cloud ice'
1607 end if
1608 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
1609 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1610 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1611 atmosphere(1)%cloud(2)%effective_radius(k) = 75.
1612 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1613 rhox=1000.
1614 rho=pmid(i,j,k)/(rd*t(i,j,k)*(1.+d608*q(i,j,k)))
1615 atmosphere(1)%cloud(3)%effective_radius(k) = 0.
1616 if(nrain(i,j,k)>0.) &
1617 atmosphere(1)%cloud(3)%effective_radius(k) = &
1618 1.0e6*1.5*(rho*qqr(i,j,k)/(pi*rhox*nrain(i,j,k)))**(1./3.)
1619 atmosphere(1)%cloud(3)%water_content(k) = max(0.,qqr(i,j,k)*dpovg)
1620 if(f_rimef(i,j,k)<=5.0)then
1621 rhox=100
1622 if(nlice(i,j,k)>0.) &
1623 atmosphere(1)%cloud(4)%effective_radius(k) = &
1624 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.) !convert to microns
1625 atmosphere(1)%cloud(4)%water_content(k) = max(0.,qqs(i,j,k)*dpovg)
1626 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1627 atmosphere(1)%cloud(5)%water_content(k) =0.
1628 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1629 atmosphere(1)%cloud(6)%water_content(k) =0.
1630 else if(f_rimef(i,j,k)<=20.0)then
1631 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1632 atmosphere(1)%cloud(4)%water_content(k) =0.
1633 rhox=400.
1634 if(nlice(i,j,k)>0.) &
1635 atmosphere(1)%cloud(5)%effective_radius(k) = &
1636 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1637 atmosphere(1)%cloud(5)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1638 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1639 atmosphere(1)%cloud(6)%water_content(k) =0.
1640 else
1641 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1642 atmosphere(1)%cloud(4)%water_content(k) =0.
1643 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1644 atmosphere(1)%cloud(5)%water_content(k) =0.
1645 rhox=900.
1646 if(nlice(i,j,k)>0.) &
1647 atmosphere(1)%cloud(6)%effective_radius(k) = &
1648 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1649 atmosphere(1)%cloud(6)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1650 end if
1651 if(debugprint .and. i==im/2 .and. j==jsta) &
1652 print*,'sample precip ice radius= ',i,j,k, f_rimef(i,j,k), &
1653 atmosphere(1)%cloud(4)%effective_radius(k), atmosphere(1)%cloud(4)%water_content(k), &
1654 atmosphere(1)%cloud(5)%effective_radius(k), atmosphere(1)%cloud(5)%water_content(k), &
1655 atmosphere(1)%cloud(6)%effective_radius(k), atmosphere(1)%cloud(6)%water_content(k)
1656
1657 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 .or. &
1658 imp_physics==28 .or. imp_physics==11)then
1659 atmosphere(1)%cloud(1)%water_content(k)=max(0.,qqw(i,j,k)*dpovg)
1660 atmosphere(1)%cloud(2)%water_content(k)=max(0.,qqi(i,j,k)*dpovg)
1661 atmosphere(1)%cloud(3)%water_content(k)=max(0.,qqr(i,j,k)*dpovg)
1662 atmosphere(1)%cloud(4)%water_content(k)=max(0.,qqs(i,j,k)*dpovg)
1663 atmosphere(1)%cloud(5)%water_content(k)=max(0.,qqg(i,j,k)*dpovg)
1664! use the effective radii output directly by the model where possible
1665 if(effrl(i,j,k)/=spval)then
1666 atmosphere(1)%cloud(1)%effective_radius(k)=min(max(effrl(i,j,k),2.5),50.)
1667 else
1668 atmosphere(1)%cloud(1)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1669 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1670 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'C')
1671 endif
1672 if(effri(i,j,k)/=spval)then
1673 atmosphere(1)%cloud(2)%effective_radius(k)=min(max(effri(i,j,k),2.5),125.)
1674 else
1675 atmosphere(1)%cloud(2)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1676 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1677 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'I')
1678 endif
1679 if(effrs(i,j,k)/=spval)then
1680 atmosphere(1)%cloud(4)%effective_radius(k)=min(max(effrs(i,j,k),2.5),1000.)
1681 else
1682 atmosphere(1)%cloud(4)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1683 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1684 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'S')
1685 endif
1686 atmosphere(1)%cloud(3)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1687 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1688 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'R')
1689 atmosphere(1)%cloud(5)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1690 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1691 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,'G')
1692 end if
1693 end do
1694!Meng 09/2018 modify two model layer having identical pressure
1695 do k = 1,lm-1
1696 if (abs(atmosphere(1)%pressure(k)-atmosphere(1)%pressure(k+1)) &
1697 < 0.005) then
1698 atmosphere(1)%pressure(k)=atmosphere(1)%pressure(k)+0.005
1699 endif
1700 enddo
1701 !bsf - start
1702 !-- Add subgrid-scale convective clouds for WRF runs
1703 if(icu_physics==2) then
1704 lcbot=nint(hbot(i,j))
1705 lctop=nint(htop(i,j))
1706 if (lcbot-lctop > 1) then
1707 !-- q_conv = assumed grid-averaged cloud water/ice condensate from convection (Cu)
1708 ! In "params" Qconv=0.1e-3 and TRAD_ice=253.15; cnvcfr is the Cu cloud fraction
1709 ! calculated as a function of Cu rain rate (Slingo, 1987) in subroutine MDLFLD
1710 q_conv = cnvcfr(i,j)*qconv
1711 do k = lctop,lcbot
1712 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g
1713 if (t(i,j,k) < trad_ice) then
1714 atmosphere(1)%cloud(2)%water_content(k) = &
1715 atmosphere(1)%cloud(2)%water_content(k) + dpovg*q_conv
1716 else
1717 atmosphere(1)%cloud(1)%water_content(k) = &
1718 atmosphere(1)%cloud(1)%water_content(k) + dpovg*q_conv
1719 endif
1720 end do !-- do k = lctop,lcbot
1721 endif !-- if (lcbot-lctop > 1) then
1722 endif !-- if (MODELNAME == 'NMM' .OR. MODELNAME == 'NCAR') then
1723 !bsf - end
1724
1725 ! call crtm forward model
1726 error_status = crtm_forward(atmosphere,surface, &
1727 geometryinfo,channelinfo(sensorindex:sensorindex), &
1728 rtsolution)
1729 if (error_status /=0) then
1730 print*,'***ERROR*** during crtm_forward call ', &
1731 error_status
1732 do n=1,channelinfo(sensorindex)%n_channels
1733 tb(i,j,n)=spval
1734 end do
1735 else
1736 do n=1,channelinfo(sensorindex)%n_channels
1737 tb(i,j,n)=rtsolution(n,1)%brightness_temperature
1738 end do
1739 if(i==ii.and.j==jj.and.debugprint) then
1740 do n=1,channelinfo(sensorindex)%n_channels
1741 3303 format('Sample rtsolution(',i0,',',i0,') in CALRAD = ',f0.3)
1742! print 3303,n,1,rtsolution(n,1)%brightness_temperature
1743 enddo
1744 do n=1,channelinfo(sensorindex)%n_channels
1745 3304 format('Sample tb(',i0,',',i0,',',i0,') in CALRAD = ',f0.3)
1746 print 3304,i,j,n,tb(i,j,n)
1747 enddo
1748 endif
1749 ! if(tb1(i,j) < 400. ) &
1750 ! & print*,'good tb1 ',i,j,tb1(i,j),gdlat(i,j),gdlon(i,j)
1751 ! if(tb2(i,j) > 400.)print*,'bad tb2 ',i,j,tb2(i,j)
1752 ! if(tb3(i,j) > 400.)print*,'bad tb3 ',i,j,tb3(i,j)
1753 ! if(tb4(i,j) > 400.)print*,'bad tb4 ',i,j,tb4(i,j)
1754 end if
1755 else
1756 do n=1,channelinfo(sensorindex)%n_channels
1757 tb(i,j,n)=spval
1758 end do
1759 END IF ! endif block for allowable satellite zenith angle
1760 end do loopi2 ! end loop for i
1761 end do ! end loop for j
1762
1763 ! error_status = crtm_destroy(channelinfo)
1764 ! if (error_status /= success) &
1765 ! & print*,'ERROR*** crtm_destroy error_status=',error_status
1766
1767 if (isis=='ssmi_f13')then ! writing ssmi to grib (37 & 85 GHz)
1768 nc=0
1769 do ixchan=1,7
1770 ichan=ixchan
1771 igot=iget(800)
1772 if(igot>0) then
1773 if(lvls(ixchan,igot)==1)then
1774 nc=nc+1
1775 do j=jsta,jend
1776 do i=ista,iend
1777 grid1(i,j)=tb(i,j,nc)
1778 enddo
1779 enddo
1780 if (grib=="grib2") then
1781 cfld=cfld+1
1782 fld_info(cfld)%ifld=iavblfld(igot)
1783 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1784 endif
1785 endif
1786 endif
1787 enddo
1788 end if ! end of outputting ssmi f13
1789 if (isis=='ssmi_f14')then ! writing ssmi to grib (19,37 & 85 GHz)
1790 nc=0
1791 do ixchan=1,7
1792 ichan=ixchan
1793 igot=iget(806)
1794 if(igot>0) then
1795 if(lvls(ixchan,igot)==1)then
1796 nc=nc+1
1797 do j=jsta,jend
1798 do i=ista,iend
1799 grid1(i,j)=tb(i,j,nc)
1800 enddo
1801 enddo
1802 if (grib=="grib2") then
1803 cfld=cfld+1
1804 fld_info(cfld)%ifld=iavblfld(igot)
1805 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1806 endif
1807 endif
1808 endif
1809 enddo
1810 end if ! end of outputting ssmi f14
1811 if (isis=='ssmi_f15')then ! writing ssmi to grib (19,37 & 85 GHz)
1812 nc=0
1813 do ixchan=1,7
1814 ichan=ixchan
1815 igot=iget(812)
1816 if(igot>0) then
1817 if(lvls(ixchan,igot)==1)then
1818 nc=nc+1
1819 do j=jsta,jend
1820 do i=ista,iend
1821 grid1(i,j)=tb(i,j,nc)
1822 enddo
1823 enddo
1824 if (grib=="grib2") then
1825 cfld=cfld+1
1826 fld_info(cfld)%ifld=iavblfld(igot)
1827 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1828 endif
1829 endif
1830 endif
1831 enddo
1832 end if ! end of outputting ssmi f15
1833 if (isis=='ssmis_f16')then ! writing ssmis to grib (183,19,37 & 85GHz)
1834 nc=0
1835 do ixchan=1,24
1836 ichan=ixchan
1837 igot=iget(818)
1838 if(igot>0) then
1839 print*,'ixchan,lvls=',ixchan,lvls(ixchan,igot)
1840 if(lvls(ixchan,igot)==1)then
1841 nc=nc+1
1842 do j=jsta,jend
1843 do i=ista,iend
1844 grid1(i,j)=tb(i,j,nc)
1845 enddo
1846 enddo
1847 if (grib=="grib2") then
1848 cfld=cfld+1
1849 fld_info(cfld)%ifld=iavblfld(igot)
1850 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1851 endif
1852 endif
1853 endif
1854 enddo
1855 end if ! end of outputting ssmis f16
1856 if(isis=='ssmis_f17') then ! writing ssmis f17 to grib (37, 91GHz)
1857 do ixchan=1,4
1858 ichan=14+ixchan
1859 igot=iget(824+ixchan)
1860 if(igot>0)then
1861 do j=jsta,jend
1862 do i=ista,iend
1863 grid1(i,j)=tb(i,j,ichan)
1864 enddo
1865 enddo
1866 if(grib=="grib2" )then
1867 cfld=cfld+1
1868 fld_info(cfld)%ifld=iavblfld(igot)
1869 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1870 endif
1871 endif
1872 enddo
1873 endif ! end of outputting ssmis f17
1874 if (isis=='ssmis_f18')then ! writing ssmis to grib (183,19,37 &85GHz)
1875 nc=0
1876 do ixchan=1,24
1877 ichan=ixchan
1878 igot=iget(832)
1879 if(igot>0) then
1880 if(lvls(ixchan,igot)==1)then
1881 nc=nc+1
1882 do j=jsta,jend
1883 do i=ista,iend
1884 grid1(i,j)=tb(i,j,nc)
1885 enddo
1886 enddo
1887 if (grib=="grib2") then
1888 cfld=cfld+1
1889 fld_info(cfld)%ifld=iavblfld(igot)
1890 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1891 endif
1892 endif
1893 endif
1894 enddo
1895 end if ! end of outputting ssmis f18
1896 if (isis=='ssmis_f19')then ! writing ssmis to grib (183,19,37 &85GHz)
1897 nc=0
1898 do ixchan=1,24
1899 ichan=ixchan
1900 igot=iget(839)
1901 if(igot>0) then
1902 if(lvls(ixchan,igot)==1)then
1903 nc=nc+1
1904 do j=jsta,jend
1905 do i=ista,iend
1906 grid1(i,j)=tb(i,j,nc)
1907 enddo
1908 enddo
1909 if (grib=="grib2") then
1910 cfld=cfld+1
1911 fld_info(cfld)%ifld=iavblfld(igot)
1912 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1913 endif
1914 endif
1915 endif
1916 enddo
1917 end if ! end of outputting ssmis f19
1918 if (isis=='ssmis_f20')then ! writing ssmis to grib (183,19,37 &85GHz)
1919 nc=0
1920 do ixchan=1,24
1921 ichan=ixchan
1922 igot=iget(846)
1923 if(igot>0) then
1924 if(lvls(ixchan,igot)==1)then
1925 nc=nc+1
1926 do j=jsta,jend
1927 do i=ista,iend
1928 grid1(i,j)=tb(i,j,nc)
1929 enddo
1930 enddo
1931 if (grib=="grib2") then
1932 cfld=cfld+1
1933 fld_info(cfld)%ifld=iavblfld(igot)
1934 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1935 endif
1936 endif
1937 endif
1938 enddo
1939 end if ! end of outputting ssmis f20
1940 if(isis=='imgr_mt2') then ! writing MTSAT-2 to grib
1941 nc=0
1942 do ichan=1,4
1943 igot=iget(860)
1944 if(lvls(ichan,igot)==1)then
1945 nc=nc+1
1946 do j=jsta,jend
1947 do i=ista,iend
1948 grid1(i,j)=tb(i,j,nc)
1949 enddo
1950 enddo
1951 if(grib=="grib2") then
1952 cfld=cfld+1
1953 fld_info(cfld)%ifld=iavblfld(igot)
1954 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1955 endif
1956 endif
1957 enddo
1958 endif
1959 if(isis=='imgr_mt1r') then ! writing MTSAT-1r to grib
1960 nc=0
1961 do ichan=1,4
1962 igot=iget(864)
1963 if(lvls(ichan,igot)==1)then
1964 nc=nc+1
1965 do j=jsta,jend
1966 do i=ista,iend
1967 grid1(i,j)=tb(i,j,nc)
1968 enddo
1969 enddo
1970 if(grib=="grib2" )then
1971 cfld=cfld+1
1972 fld_info(cfld)%ifld=iavblfld(igot)
1973 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1974 endif
1975 endif
1976 enddo
1977 endif
1978 if_insat3d: if(isis=='imgr_insat3d') then ! writing MTSAT-1r to grib
1979 nc=0
1980 do ichan=1,4
1981 igot=iget(865)
1982 if(lvls(ichan,igot)==1)then
1983 nc=nc+1
1984 do j=jsta,jend
1985 do i=ista,iend
1986 grid1(i,j)=tb(i,j,nc)
1987 enddo
1988 enddo
1989 if(grib=="grib2" )then
1990 cfld=cfld+1
1991 fld_info(cfld)%ifld=iavblfld(igot)
1992 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1993 endif
1994 endif
1995 enddo
1996 endif if_insat3d
1997 if (isis=='imgr_g11')then ! writing goes 11 to grib
1998 do ixchan=1,4
1999 ichan=ixchan
2000 igot=iget(459+ixchan)
2001 if(igot>0) then
2002 do j=jsta,jend
2003 do i=ista,iend
2004 grid1(i,j)=tb(i,j,ichan)
2005 enddo
2006 enddo
2007 if(grib=="grib2" )then
2008 cfld=cfld+1
2009 fld_info(cfld)%ifld=iavblfld(igot)
2010 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2011 endif
2012 endif
2013 enddo
2014 endif ! end of outputting goes 11
2015 if (isis=='imgr_g12')then ! writing goes 12 to grib
2016 do ixchan=1,4
2017 ichan=ixchan
2018 igot=iget(455+ixchan)
2019 if(igot>0) then
2020 do j=jsta,jend
2021 do i=ista,iend
2022 grid1(i,j)=tb(i,j,ichan)
2023 enddo
2024 enddo
2025 if(grib=="grib2" )then
2026 cfld=cfld+1
2027 fld_info(cfld)%ifld=iavblfld(igot)
2028 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2029 endif
2030 endif
2031 enddo
2032 end if ! end of outputting goes 12
2033 if (isis=='seviri_m10')then ! writing msg/severi 10
2034 nc=0
2035 do ixchan=1,8
2036 ichan=ixchan
2037 igot=iget(876)
2038 if(igot>0) then
2039 if(lvls(ixchan,igot)==1)then
2040 nc=nc+1
2041 do j=jsta,jend
2042 do i=ista,iend
2043 grid1(i,j)=tb(i,j,nc)
2044 enddo
2045 enddo
2046 if (grib=="grib2") then
2047 cfld=cfld+1
2048 fld_info(cfld)%ifld=iavblfld(igot)
2049 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2050 endif
2051 endif
2052 endif
2053 enddo
2054 end if ! end of outputting msg/seviri 10
2055 if (isis=='imgr_g13')then ! writing goes 13 to grib
2056 nc=0
2057 do ixchan=1,4
2058 ichan=ixchan
2059 igot=iget(868)
2060 if(igot>0) then
2061 if(lvls(ixchan,igot)==1)then
2062 nc=nc+1
2063 do j=jsta,jend
2064 do i=ista,iend
2065 grid1(i,j)=tb(i,j,nc)
2066 enddo
2067 enddo
2068 if (grib=="grib2") then
2069 cfld=cfld+1
2070 fld_info(cfld)%ifld=iavblfld(igot)
2071 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2072 endif
2073 endif
2074 endif
2075 enddo
2076 end if ! end of outputting goes 13
2077 if (isis=='imgr_g15')then ! writing goes 15 to grib
2078 nc=0
2079 do ixchan=1,4
2080 ichan=ixchan
2081 igot=iget(872)
2082 if(igot>0) then
2083 if(lvls(ixchan,igot)==1)then
2084 nc=nc+1
2085 do j=jsta,jend
2086 do i=ista,iend
2087 grid1(i,j)=tb(i,j,nc)
2088 enddo
2089 enddo
2090 if (grib=="grib2") then
2091 cfld=cfld+1
2092 fld_info(cfld)%ifld=iavblfld(igot)
2093 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2094 endif
2095 endif
2096 endif
2097 enddo
2098 end if ! end of outputting goes 15
2099 if (isis=='abi_g16')then ! writing goes 16 to grib
2100 nc=0
2101 do ixchan=1,10
2102 ichan=ixchan
2103 igot=iget(926+ixchan)
2104 if(igot>0)then
2105 do j=jsta,jend
2106 do i=ista,iend
2107 grid1(i,j)=tb(i,j,ichan)
2108 enddo
2109 enddo
2110 if(grib=="grib2" )then
2111 cfld=cfld+1
2112 fld_info(cfld)%ifld=iavblfld(igot)
2113 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2114 endif
2115 endif
2116 enddo ! channel loop
2117 end if ! end of outputting goes 16
2118 if (isis=='abi_g17')then ! writing goes 17 to grib
2119 nc=0
2120 do ixchan=1,10
2121 ichan=ixchan
2122 igot=iget(936+ixchan)
2123 if(igot>0)then
2124 do j=jsta,jend
2125 do i=ista,iend
2126 grid1(i,j)=tb(i,j,ichan)
2127 enddo
2128 enddo
2129 if(grib=="grib2" )then
2130 cfld=cfld+1
2131 fld_info(cfld)%ifld=iavblfld(igot)
2132 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2133 endif
2134 endif
2135 enddo ! channel loop
2136 end if ! end of outputting goes 17
2137! Wm Lewis updated idx for g18 on 3 JUN 2022
2138 if (isis=='abi_g18')then ! writing goes 18 to grib
2139 nc=0
2140 do ixchan=1,10
2141 igot=iget(530+ixchan)
2142 ichan=ixchan
2143 if(igot>0)then
2144 do j=jsta,jend
2145 do i=1,im
2146 grid1(i,j)=tb(i,j,ichan)
2147 enddo
2148 enddo
2149 id(1:25) = 0
2150 id(02) = 2
2151 id(08) = 118
2152 id(09) = 109
2153 if(grib=="grib2" )then
2154 cfld=cfld+1
2155 fld_info(cfld)%ifld=iavblfld(igot)
2156 datapd(1:im,1:jend-jsta+1,cfld)=grid1(1:im,jsta:jend)
2157 endif
2158 endif
2159 enddo ! channel loop
2160 end if ! end of outputting goes 18
2161 if(isis=='ahi_himawari8') then ! writing Himawari-8 AHI to grib
2162 do ichan=1,10
2163 igot=iget(968+ichan)
2164 if(igot>0)then
2165 do j=jsta,jend
2166 do i=ista,iend
2167 grid1(i,j)=tb(i,j,ichan)
2168 enddo
2169 enddo
2170 if(grib=="grib2" )then
2171 cfld=cfld+1
2172 fld_info(cfld)%ifld=iavblfld(igot)
2173 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2174 endif
2175 endif
2176 enddo
2177 endif ! end of outputting himawari-8 ahi
2178
2179 end if nonnadir ! end if for computing simulated radiance with zenith angle correction
2180
2181 ! Deallocate arrays
2182 CALL crtm_atmosphere_destroy(atmosphere(1))
2183 CALL crtm_surface_destroy(surface(1))
2184 CALL crtm_rtsolution_destroy(rtsolution)
2185 if (crtm_atmosphere_associated(atmosphere(1))) &
2186 write(6,*)' ***ERROR** destroying atmosphere.'
2187 if (crtm_surface_associated(surface(1))) &
2188 write(6,*)' ***ERROR** destroying surface.'
2189 if (any(crtm_rtsolution_associated(rtsolution))) &
2190 write(6,*)' ***ERROR** destroying rtsolution.'
2191 deallocate(tb)
2192 deallocate (rtsolution)
2193
2194!
2195 end if sensor_avail ! end of if block for only calling crtm when sepcific sensor is requested in the control file
2196 end do sensordo ! end looping for different satellite
2197 error_status = crtm_destroy(channelinfo)
2198 if (error_status /= success) &
2199 print*,'ERROR*** crtm_destroy error_status=',error_status
2200 deallocate(channelinfo)
2201 if (allocated(model_to_crtm)) deallocate(model_to_crtm)
2202 endif ifactive ! for all iget logical
2203 return
2204end SUBROUTINE calrad_wcloud
2205
2206!-------------------------------------------------------------------------------
2226!-------------------------------------------------------------------------------
2227
2228REAL FUNCTION effr(pmid,t,q,qqw,qqi,qqr,f_rimef, nlice, nrain, &
2229 qqs,qqg,qqnr,qqni,qqnw,mp_opt,species)
2230
2231! JASON OTKIN AND WILLIAM LEWIS
2232! 09 DECEMBER 2014
2233! Greg Thompson, 20200924
2234
2235 use params_mod, only: pi, rd, d608, rg
2236
2237 implicit none
2238
2239 real :: pmid,t,q,qqw,qqi,qqr,qqs,qqg,f_rimef,nlice,nrain
2240 real :: qqnr,qqni,qqnw
2241 character(LEN=1) :: species
2242
2243 integer :: n,count,count1,mp_opt
2244 real :: rho, ncc, rhox
2245 real :: n0_s, n0_r, n0_g
2246 real :: lambdar, lambdas, lambdag
2247
2248!-------------------------------------------------------------------------------
2249! GAMMA FUNCTION & RELATED VARIABLES
2250!-------------------------------------------------------------------------------
2251
2252 real :: gamma
2253 real :: gamma_crg, gamma_s
2254! real :: gamma_i
2255
2256 real :: wgamma, gammln
2257
2258 real :: rc,am_c,bm_c,cce(3,15),ccg(3,15),ocg1(15),ocg2(15)
2259 integer :: nu_c
2260
2261 real, dimension(0:15), parameter:: g_ratio = (/6,24,60,120,210, &
2262 & 336,504,720,990,1320,1716,2184,2730,3360,4080,4896/)
2263
2264 real :: rr, mu_r, am_r, bm_r, cre(3), crg(3), ore1, org1, org2
2265 real :: mvd_r, ron_sl, ron_r0, ron_c0, ron_c1, ron_c2, obmr
2266
2267 real :: ri, mu_i, am_i, bm_i, cie(3), cig(3), oig1, oig2, obmi
2268
2269 real :: rs, am_s, oams, cse(3)
2270 real :: loga, a, b, tc0, smob, smo2, smoc
2271 REAL, PARAMETER:: mu_s = 0.6357
2272 REAL, PARAMETER:: kap0 = 490.6
2273 REAL, PARAMETER:: kap1 = 17.46
2274 REAL, PARAMETER:: lam0 = 20.78
2275 REAL, PARAMETER:: lam1 = 3.29
2276
2277!-------------------------------------------------------------------------------
2278! MINIMUM/MAXIMUM CONSTANTS FOR ALL SCHEMES
2279!-------------------------------------------------------------------------------
2280
2281 real, parameter :: eps=0.622, beta_crg=3., beta_i=2.,beta_s=2.4
2282
2283 real, parameter :: min_qc=1.e-7, min_qr=1.e-7, min_qi=1.e-8,min_qs=1.e-8, min_qg=1.e-7
2284 real, parameter :: min_c=2.e-6, min_r=20.e-6, min_i=4.e-6,min_s=20.e-6, min_g=20.e-6
2285 real, parameter :: max_c=1.e-2, max_r=1.e-2, max_i=1.e-3,max_s=2.e-2, max_g=5.e-0
2286
2287 real :: am_g, bm_g, mu_g
2288 real :: cgg(3), cge(3), oge1, obmg, ogg1, ogg2
2289
2290 real :: ygra1, zans1, rg2
2291 double precision :: no_exp, no_min, lm_exp, lamg, lamc, lamr, lami, lams
2292
2293!-------------------------------------------------------------------------------
2294! WSM6-SPECIFIC ARRAYS
2295!-------------------------------------------------------------------------------
2296
2297 real :: wsm6_nci, xmi, xmitemp
2298
2299!-------------------------------------------------------------------------------
2300! CONSTANTS FOR WSM6 MICROPHYSICS SCHEME
2301!-------------------------------------------------------------------------------
2302
2303 real, parameter :: wsm6_cnp=3.e8, wsm6_rhor=1000.
2304 real, parameter :: wsm6_rhos=100., wsm6_rhog=500.
2305 real, parameter :: wsm6_dimax=500.e-6, wsm6_dicon=11.9
2306 real, parameter :: wsm6_alpha=.12, wsm6_t0c=273.15
2307 real, parameter :: wsm6_n0s=2.e6, wsm6_n0smax=1.e11
2308 real, parameter :: wsm6_n0r=8.e6, wsm6_n0g=4.e6
2309 real, parameter :: wsm6_qmin=1.e-15
2310
2311!-------------------------------------------------------------------------------
2312! CONSTANTS FOR LIN MICROPHYSICS SCHEME
2313!-------------------------------------------------------------------------------
2314
2315 real, parameter :: lin_rhoi=100., lin_rhor=1000., lin_rhos=100.
2316 real, parameter :: lin_rhog=400., lin_cnp=3.e8
2317 real, parameter :: lin_n0r=8.e6, lin_n0s=3.e6, lin_n0g=4.e6
2318
2319!-------------------------------------------------------------------------------
2320! CONSTANTS FOR NEW THOMPSON MICROPHYSICS SCHEME (FOR WRF VERSIONS 3.1 AND UP)
2321!-------------------------------------------------------------------------------
2322
2323 real, parameter :: nthom_rhor=1000., nthom_rhos=100.
2324! WM LEWIS updated rhog to 500 from 400
2325 real, parameter :: nthom_rhog=500., nthom_rhoi=890.
2326 real, parameter :: nthom_gon_min=1.e4, nthom_gon_max=3.e6
2327 real, parameter :: nthom_nt_c=100.e6
2328
2329 real, parameter :: nthom_min_nci=5.e2
2330 real, parameter :: nthom_min_ncr=1.e-6
2331
2332 real, parameter :: nthom_bm_s=2.0 !this is important
2333
2334 real :: nci2, ncc2, ncr2
2335
2336 real, dimension(10), parameter :: &
2337 nthom_sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
2338 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
2339 real, dimension(10), parameter :: &
2340 nthom_sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
2341 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
2342
2343!-------------------------------------------------------------------------------
2344! CONSTANTS FOR GFDL MICROPHYSICS SCHEME - which is Lin for precip clouds
2345!-------------------------------------------------------------------------------
2346
2347 real, parameter :: gfdl_rhoi=100., gfdl_rhor=1000., gfdl_rhos=100.
2348 real, parameter :: gfdl_rhog=400., gfdl_cnp=3.e8
2349 real, parameter :: gfdl_tice = 273.16
2350
2351 real, parameter :: gfdl_qmin = 1.0e-5, gfdl_ccn = 1.0e8, gfdl_beta = 1.22
2352 real, parameter :: gfdl_gammar = 17.837789, gfdl_gammas = 8.2850630, gfdl_gammag = 11.631769
2353 real, parameter :: gfdl_alphar = 0.8, gfdl_alphas = 0.25, gfdl_alphag = 0.5
2354 real, parameter :: gfdl_n0r=8.e6, gfdl_n0s=3.e6, gfdl_n0g=4.e6
2355
2356 real, parameter :: gfdl_rewmin = 5.0, gfdl_rewmax = 10.0
2357 real, parameter :: gfdl_reimin = 10.0, gfdl_reimax = 150.0
2358 real, parameter :: gfdl_rermin = 0.0, gfdl_rermax = 10000.0
2359 real, parameter :: gfdl_resmin = 0.0, gfdl_resmax = 10000.0
2360 real, parameter :: gfdl_regmin = 0.0, gfdl_regmax = 10000.0
2361
2362
2363
2364 if(mp_opt==6) then !WSM6 SCHEME
2365
2366 n0_r = wsm6_n0r
2367 n0_g = wsm6_n0g
2368 n0_s = wsm6_n0s
2369
2370 elseif(mp_opt==2)then !LIN SCHEME
2371
2372 n0_r = lin_n0r
2373 n0_g = lin_n0g
2374 n0_s = lin_n0s
2375
2376 endif
2377
2378 gamma_crg = 6.0 ! gamma(1.0 + beta_crg)
2379 gamma_s = 2.981134 ! gamma(1.0 + beta_s)
2380! gamma_i = 2.0 ! gamma(1.0 + beta_i)
2381
2382!------------------------------------------------------------------------------
2383! SET DIAMETER ARRAYS TO ZERO, COMPUTE DENSITY
2384!------------------------------------------------------------------------------
2385
2386 effr=0.
2387
2388 rho=pmid/(rd*t*(1.+d608*q))
2389
2390
2391 if(mp_opt==6)then
2392
2393 SELECT CASE(species)
2394
2395 CASE("C")
2396
2397 if ( qqw>min_qc ) then !cloud diameter: assume constant # concentration
2398 effr = 1.0e6*(( 6. * rho * qqw ) / &
2399 (pi * wsm6_rhor * wsm6_cnp))**(1/3.)
2400
2401 endif
2402
2403 CASE("R")
2404
2405 if ( qqr>min_qr ) then !rain diameter: assume gamma distribution
2406 effr = 1.0e6*( ( 6. * rho * qqr ) / &
2407 ( pi * wsm6_rhor * n0_r * gamma_crg ) ) ** (1/(1+beta_crg ) )
2408 endif
2409
2410 CASE("G")
2411
2412 if ( qqg>min_qg ) then !graupel diameter: assume gamma distribution
2413 effr = 1.0e6*( ( 6. * rho * qqg ) / &
2414 ( pi * wsm6_rhog * n0_g * gamma_crg ) ) ** (1/(1+beta_crg ) )
2415 endif
2416
2417 CASE("S")
2418
2419 if ( qqs>min_qs ) then !snow diameter: assume gamma distribution
2420 effr = 1.0e6*( ( 6. * rho * qqs ) / &
2421 ( pi * wsm6_rhos * n0_s * gamma_s ) ) ** ( 1/(1+beta_s) )
2422 endif
2423
2424! ICE DIAMETER: CALCULATED USING METHOD OUTLINED IN WRF BROWSER. Refer to
2425! phys/module_mp_wsm6.F (Vice:fallout of ice crystal).
2426
2427 CASE("I")
2428
2429 if ( qqi>min_qi ) then !ice diameter
2430! wsm6_nci = min(max(5.38e7*(rho*max(qqi,wsm6_qmin)),1.e3),1.e6)
2431! xmi = rho * qqi / wsm6_nci
2432! effr = 1.0E6*min( sqrt(xmi), wsm6_dimax)
2433!! from wsm6, HWRF ver 3.6:
2434!! temp = (den(i,k)*max(qci(i,k,2),qmin))
2435!! temp = sqrt(sqrt(temp*temp*temp))
2436!! xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2437!! diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
2438 xmitemp=rho*max(qqi,wsm6_qmin)
2439 xmitemp=sqrt(sqrt(xmitemp*xmitemp*xmitemp))
2440 xmi= min(max(5.38e7*xmitemp,1.e3),1.e6)
2441 effr = 1.0e6*max(min(wsm6_dicon * sqrt(xmi),wsm6_dimax), 1.e-25)
2442 endif
2443
2444 END SELECT
2445
2446 elseif(mp_opt==2)then
2447
2448 SELECT CASE(species)
2449
2450 CASE("C")
2451
2452 if ( qqw > min_qc ) then !cloud diameter: assume constant # concentration
2453 effr = 1.0e6*(( 6. * rho * qqw ) / &
2454 (pi * lin_rhor * lin_cnp))**(1/3.)
2455 endif
2456
2457 CASE("R")
2458
2459 if ( qqr > min_qr ) then !rain diameter: assume gamma distribution
2460 effr = 1.0e6*( ( 6. * rho * qqr ) / &
2461 ( pi * lin_rhor * n0_r * gamma_crg ) ) ** (1/(1+beta_crg ) )
2462 endif
2463
2464 CASE("I")
2465
2466 if ( qqi > min_qi ) then !ice diameter: assume constant # concentrtion
2467 effr = 1.0e6*( ( 6. * rho * qqi ) / &
2468 ( pi * lin_rhoi * lin_cnp ) ) ** ( 1/3.)
2469 endif
2470
2471 CASE("S")
2472
2473 if ( qqs > min_qs ) then !snow diameter: assume gamma distribution
2474 effr = 1.0e6*( ( 6. * rho * qqs ) / &
2475 ( pi * lin_rhos * n0_s * gamma_s ) ) ** ( 1/(1+beta_s) )
2476 endif
2477
2478 CASE("G")
2479
2480 if ( qqg > min_qg ) then !graupel diameter: assume gamma distribution
2481 effr = 1.0e6*( ( 6. * rho * qqg ) / &
2482 ( pi * lin_rhog * n0_g * gamma_crg ) ) ** (1/(1+beta_crg ) )
2483 endif
2484
2485 END SELECT
2486
2487 elseif(mp_opt==8 .or. mp_opt==28)then
2488
2489! rain section
2490
2491 bm_r = 3.0
2492 mu_r = 0.0
2493 obmr = 1.0 / bm_r
2494 am_r = pi * nthom_rhor / 6.0
2495
2496 cre(1) = bm_r + 1.
2497 cre(2) = mu_r + 1.
2498 cre(3) = bm_r + mu_r + 1.
2499
2500 crg(1) = wgamma(cre(1))
2501 crg(2) = wgamma(cre(2))
2502 crg(3) = wgamma(cre(3))
2503
2504 ore1 = 1. / cre(1)
2505 org1 = 1. / crg(1)
2506 org2 = 1. / crg(2)
2507
2508! cloud section
2509
2510 bm_c = bm_r
2511
2512 do n = 1, 15
2513 cce(1,n) = n + 1. ! Substitute variable value of mu_c
2514 cce(2,n) = bm_c + n + 1. ! Substitute variable value of mu_c
2515
2516 ccg(1,n) = wgamma(cce(1,n))
2517 ccg(2,n) = wgamma(cce(2,n))
2518
2519 ocg1(n) = 1./ccg(1,n)
2520 ocg2(n) = 1./ccg(2,n)
2521 enddo
2522
2523! ice section
2524
2525 am_i = pi * nthom_rhoi / 6.0
2526 bm_i = 3.0
2527 mu_i = 0.
2528
2529 cie(1) = mu_i + 1.
2530 cie(2) = bm_i + mu_i + 1.
2531
2532 cig(1) = wgamma(cie(1))
2533 cig(2) = wgamma(cie(2))
2534
2535 oig1 = 1./cig(1)
2536 oig2 = 1./cig(2)
2537 obmi = 1./bm_i
2538
2539! snow section
2540
2541 am_s = 0.069
2542
2543 oams = 1./am_s
2544
2545 cse(1) = nthom_bm_s + 1.
2546
2547! graupel section
2548 bm_g = 3.0
2549 mu_g = 0.0
2550 obmg = 1.0 / bm_g
2551 am_g = pi * nthom_rhog / 6.0
2552
2553 cge(1) = bm_g + 1.
2554 cge(2) = mu_g + 1.
2555 cge(3) = bm_g + mu_g + 1.
2556
2557 cgg(1) = wgamma(cge(1))
2558 cgg(2) = wgamma(cge(2))
2559 cgg(3) = wgamma(cge(3))
2560
2561 oge1 = 1. / cge(1)
2562 ogg1 = 1. / cgg(1)
2563 ogg2 = 1. / cgg(2)
2564
2565!CLOUD DIAMETER CALCULATIONS
2566
2567 SELECT CASE (species)
2568
2569 CASE("C")
2570
2571 if(qqw >= min_qc) then
2572
2573 rc = max(1.e-12, qqw * rho)
2574
2575 if (mp_opt==8) then
2576 ncc2 = nthom_nt_c
2577 elseif (mp_opt==28) then
2578 ncc2 = max(1.e-6, qqnw * rho)
2579 endif
2580
2581 if (ncc2 < 10.e6) then
2582 nu_c = 15
2583 else
2584 nu_c = min(15, nint(1000.e6/ncc2) + 2)
2585 endif
2586
2587 lamc = (ncc2/rc)**obmr * (am_r*g_ratio(nu_c))**obmr
2588
2589 effr = 1.0e6*max(4.01e-6, min(sngl(1.0d0*dble(3.+nu_c)/lamc),50.e-6))
2590
2591! old UPP
2592! effr = 2.*10.
2593
2594 endif
2595
2596!RAIN DIAMETER CALCULATIONS
2597
2598 CASE("R")
2599
2600 if( qqr > min_qr) then
2601
2602 rr = max(1.e-12, qqr * rho)
2603 ncr2 = max(1.e-6, qqnr * rho)
2604 lamr = (ncr2/rr)**obmr * (am_r*crg(3)*org2)**obmr
2605
2606 effr = 1.0e6*max(50.01e-6, min(sngl(1.0d0*dble(3.+mu_r)/lamr),1999.e-6))
2607
2608! old UPP
2609! effr=2.*200.
2610
2611! print*,'effr_rain=',effr/2.
2612
2613 endif
2614
2615!ICE DIAMETER CACLULATIONS
2616
2617 CASE("I")
2618
2619 if(qqi >= min_qi) then
2620
2621 ri = max(1.e-12, qqi * rho)
2622 nci2 = max(1.e-6, qqni * rho)
2623
2624 lami = (nci2/ri)**obmi * (am_i*cig(2)*oig1)**obmi
2625
2626 effr = 1.0e6*max(10.01e-6, min(sngl(1.0d0*dble(3.+mu_i)/lami),250.e-6))
2627
2628! old UPP
2629! effr=2.*25.
2630
2631 endif
2632
2633!SNOW DIAMETER CALCULATIONS
2634
2635 CASE("S")
2636
2637 rs = max(1.e-12, qqs * rho)
2638
2639 if(qqs >= min_qs) then
2640
2641 tc0 = min(-0.1, t-273.15)
2642 smob = rs*oams
2643
2644 if (nthom_bm_s>(2.0-1.e-3) .and. nthom_bm_s<(2.0+1.e-3))then
2645 smo2 = smob
2646 else
2647 loga = nthom_sa(1) + nthom_sa(2)*tc0 + nthom_sa(3)*nthom_bm_s+ &
2648 nthom_sa(4)*tc0*nthom_bm_s + nthom_sa(5)*tc0*tc0 +&
2649 nthom_sa(6)*nthom_bm_s*nthom_bm_s +nthom_sa(7)*tc0*tc0*nthom_bm_s + &
2650 nthom_sa(8)*tc0*nthom_bm_s*nthom_bm_s +nthom_sa(9)*tc0*tc0*tc0 + &
2651 nthom_sa(10)*nthom_bm_s*nthom_bm_s*nthom_bm_s
2652
2653 a = 10.0**loga
2654
2655 b = nthom_sb(1) + nthom_sb(2)*tc0 + nthom_sb(3)*nthom_bm_s+ &
2656 nthom_sb(4)*tc0*nthom_bm_s + nthom_sb(5)*tc0*tc0 +&
2657 nthom_sb(6)*nthom_bm_s*nthom_bm_s +nthom_sb(7)*tc0*tc0*nthom_bm_s + &
2658 nthom_sb(8)*tc0*nthom_bm_s*nthom_bm_s +nthom_sb(9)*tc0*tc0*tc0 + &
2659 nthom_sb(10)*nthom_bm_s*nthom_bm_s*nthom_bm_s
2660 smo2 = (smob/a)**(1./b)
2661 endif
2662
2663 !Calculate bm_s+1 (th) moment. Useful for diameter calcs.
2664 loga = nthom_sa(1) + nthom_sa(2)*tc0 + nthom_sa(3)*cse(1) +&
2665 nthom_sa(4)*tc0*cse(1) + nthom_sa(5)*tc0*tc0 +&
2666 nthom_sa(6)*cse(1)*cse(1) + nthom_sa(7)*tc0*tc0*cse(1)+ &
2667 nthom_sa(8)*tc0*cse(1)*cse(1) +nthom_sa(9)*tc0*tc0*tc0 + &
2668 nthom_sa(10)*cse(1)*cse(1)*cse(1)
2669
2670 a = 10.0**loga
2671
2672 b = nthom_sb(1)+ nthom_sb(2)*tc0 + nthom_sb(3)*cse(1) +&
2673 nthom_sb(4)*tc0*cse(1) + nthom_sb(5)*tc0*tc0 +&
2674 nthom_sb(6)*cse(1)*cse(1) + nthom_sb(7)*tc0*tc0*cse(1) +&
2675 nthom_sb(8)*tc0*cse(1)*cse(1) + nthom_sb(9)*tc0*tc0*tc0 + &
2676 nthom_sb(10)*cse(1)*cse(1)*cse(1)
2677
2678 smoc = a * smo2**b
2679
2680 effr = 1.0e6*max(50.e-6, min(smoc/smob, 1999.e-6))
2681
2682! print*,'snow effr=',effr
2683
2684! changing snow effr recovers "old" UPP Thompson almost exactly;
2685! i.e. the snow effr is the source of the cold discprepancy.
2686
2687! old UPP
2688! effr=2.*250.
2689
2690 endif
2691
2692 CASE("G")
2693
2694 if(qqg >= min_qg) then
2695
2696 rg2 = max(1.e-12, qqg * rho)
2697
2698 ygra1 = alog10(max(1.e-9, rg2))
2699
2700 zans1 = 3. + 2./7. * (ygra1+7.)
2701 zans1 = max(2., min(zans1, 7.))
2702
2703 no_exp = 10.**(zans1)
2704
2705 lm_exp = (no_exp*am_g*cgg(1)/rg2)**oge1
2706
2707 lamg = lm_exp * (cgg(3)*ogg2*ogg1)**obmg
2708
2709 effr= 1.0e6*max(99.e-6, min(sngl((3.0+mu_g)/lamg), 9999.e-6))
2710
2711! old UPP
2712! effr=350.
2713
2714 endif
2715
2716 END SELECT
2717
2718 elseif(mp_opt==11)then ! GFDL
2719
2720 SELECT CASE(species)
2721
2722 CASE("C")
2723
2724! cloud water (martin et al., 1994)
2725 if (qqw > min_qc) then
2726 effr = exp(1.0 / 3.0 * log((3. * qqw ) / (4. * pi * gfdl_rhor * gfdl_ccn))) * 1.0e6
2727 effr = max(gfdl_rewmin, min(gfdl_rewmax, effr))
2728 effr = effr*2. ! because need diameter here, converted to radius at exit
2729 end if
2730
2731 CASE("I")
2732
2733! cloud ice (heymsfield and mcfarquhar, 1996)
2734 if (qqi > min_qi) then
2735 if ((t-gfdl_tice) < - 50) then
2736 effr = gfdl_beta / 9.917 * exp((1 - 0.891) * log(1.0e3 * qqi)) * 1.0e3
2737 elseif ((t-gfdl_tice) < - 40.) then
2738 effr = gfdl_beta / 9.337 * exp((1 - 0.920) * log(1.0e3 * qqi)) * 1.0e3
2739 elseif ((t-gfdl_tice) < - 30.) then
2740 effr = gfdl_beta / 9.208 * exp((1 - 0.945) * log(1.0e3 * qqi)) * 1.0e3
2741 else
2742 effr = gfdl_beta / 9.387 * exp((1 - 0.969) * log(1.0e3 * qqi)) * 1.0e3
2743 endif
2744 effr = max(gfdl_reimin, min(gfdl_reimax, effr))
2745 effr = effr*2. ! because need diameter here, converted to radius at exit
2746 end if
2747
2748 CASE("R")
2749
2750 if ( qqr > min_qr ) then !rain diameter: assume gamma distribution
2751 lambdar = exp(0.25 * log(pi * gfdl_rhor * gfdl_n0r / qqr))
2752 effr = 0.5*exp(log(gfdl_gammar / 6.) / gfdl_alphar) / lambdar * 1.0e6
2753 effr = max(gfdl_rermin, min(gfdl_rermax, effr))
2754 effr = effr*2. ! because need diameter here, converted to radius at exit
2755 endif
2756
2757
2758 CASE("S")
2759
2760 if ( qqs > min_qs ) then !snow diameter: assume gamma distribution
2761 lambdas = exp(0.25 * log(pi * gfdl_rhos * gfdl_n0s / qqs))
2762 effr = 0.5 * exp(log(gfdl_gammas / 6.) / gfdl_alphas) / lambdas * 1.0e6
2763 effr = max(gfdl_resmin, min(gfdl_resmax, effr))
2764 effr = effr*2. ! because need diameter here, converted to radius at exit
2765 endif
2766
2767 CASE("G")
2768
2769 if ( qqg > min_qg ) then !graupel diameter: assume gamma distribution
2770 lambdag = exp(0.25 * log(pi * gfdl_rhog * gfdl_n0g / qqg))
2771 effr = 0.5 * exp(log(gfdl_gammag / 6.) / gfdl_alphag) / lambdag * 1.0e6
2772 effr = max(gfdl_regmin, min(gfdl_regmax, effr))
2773 effr = effr*2. ! because need diameter here, converted to radius at exit
2774 endif
2775
2776 END SELECT
2777
2778
2779 elseif(mp_opt==5.or.mp_opt==85.or.mp_opt==95)then
2780
2781 SELECT CASE (species)
2782
2783 CASE("C")
2784
2785 effr=2.*10.
2786
2787 CASE("I")
2788
2789 effr=2.*25.
2790
2791 CASE("R")
2792
2793 if( qqr > min_qr) then
2794 rhox=1000.
2795 effr=2.*1.0e6*1.5*(rho*qqr/(pi*rhox*nrain))**(1./3.)
2796
2797! old UPP
2798! effr=2.*200.
2799! effr=min(200.,effr)
2800! print*,'effr_rain=',effr/2.
2801 endif
2802
2803 CASE("S")
2804
2805 if(f_rimef<=5.0)then
2806 rhox=100.
2807 if(nlice>0.) then
2808 effr = 2.*1.0e6*1.5*(rho*qqs/(pi*rhox*nlice))**(1./3.)
2809 endif
2810 endif
2811
2812 CASE("G")
2813
2814 if(f_rimef>5.0.and.f_rimef<=20.0)then
2815 rhox=400.
2816 if(nlice>0.) then
2817 effr = 2.*1.0e6*1.5*(rho*qqs/(pi*rhox*nlice))**(1./3.)
2818 endif
2819 endif
2820
2821 CASE("H")
2822
2823 if(f_rimef>20.0)then
2824 rhox=900.
2825 if(nlice>0.) then
2826 effr = 2.*1.0e6*1.5*(rho*qqs/(pi*rhox*nlice))**(1./3.)
2827 endif
2828 endif
2829
2830 END SELECT
2831
2832
2833 endif
2834
2835!-----------------------------------------
2836! DIAMETER -> RADIUS
2837!-----------------------------------------
2838
2839 effr = 0.5*effr
2840
2841
2842end function effr
2843
2844!-------------------------------------------------------------------------------
2849!-------------------------------------------------------------------------------
2850
2851 REAL function gammln(xx)
2852! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
2853 IMPLICIT NONE
2854 REAL, INTENT(IN):: xx
2855 DOUBLE PRECISION, PARAMETER:: stp = 2.5066282746310005d0
2856 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
2857 cof = (/76.18009172947146d0, -86.50532032941677d0, &
2858 24.01409824083091d0, -1.231739572450155d0, &
2859 .1208650973866179d-2, -.5395239384953d-5/)
2860 DOUBLE PRECISION:: ser,tmp,x,y
2861 INTEGER:: j
2862
2863 x=xx
2864 y=x
2865 tmp=x+5.5d0
2866 tmp=(x+0.5d0)*log(tmp)-tmp
2867 ser=1.000000000190015d0
2868 DO 11 j=1,6
2869 y=y+1.d0
2870 ser=ser+cof(j)/y
287111 CONTINUE
2872 gammln=tmp+log(stp*ser/x)
2873 END FUNCTION gammln
2874
2875!-------------------------------------------------------------------------------
2880!-------------------------------------------------------------------------------
2881
2882 REAL function wgamma(y)
2883
2884 IMPLICIT NONE
2885 REAL, INTENT(IN):: y
2886
2887 real :: gammln
2888
2889 wgamma = exp(gammln(y))
2890
2891 END FUNCTION wgamma
2892
gammln
real function gammln(xx)
GAMMLN.
Definition CALRAD_WCLOUD_newcrtm.f:2852
calrad_wcloud
subroutine calrad_wcloud
CALRAD_WCLOUD Subroutine that computes model derived brightness temperature.
Definition CALRAD_WCLOUD_newcrtm.f:29
wgamma
real function wgamma(y)
WGAMMA.
Definition CALRAD_WCLOUD_newcrtm.f:2883
effr
real function effr(pmid, t, q, qqw, qqi, qqr, f_rimef, nlice, nrain, qqs, qqg, qqnr, qqni, qqnw, mp_opt, species)
EFFR Computes effective particle radii channel selection using LVLS from WRF_CNTRL....
Definition CALRAD_WCLOUD_newcrtm.f:2230
initpost
subroutine initpost
This routine initializes constants and variables at the start of an ETA model or post processor run.
Definition INITPOST.F:26
process
subroutine process(kth, kpv, th, pv, iostatusd3d)
process() is a driver for major post routines.
Definition PROCESS.f:40
select_channels_l
subroutine select_channels_l(channelinfo, nchannels, channels, l, igot)
select_channels_L() verifies channel information.
Definition SELECT_CHANNELS.f:70
zensun
subroutine zensun(day, time, lat, lon, pi, sun_zenith, sun_azimuth)
Makes sun zenith and sun azimuth angle.
Definition ZENSUN.f:62
ctlblk_mod::list
integer list
No longer used/supported.
Definition CTLBLK.f:104