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