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