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