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