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