30 use vrbls3d,
only: o3, pint, pmid, t, q, qqw, qqi, qqr, f_rimef, nlice, nrain, qqs, qqg, &
31 qqnr, qqni, qqnw, cfr, effri, effrl, effrs
32 use vrbls2d,
only: czen, ivgtyp, sno, pctsno, ths, vegfrc, si, u10h, v10h, u10,&
33 v10, smstot, hbot, htop, cnvcfr
34 use masks,
only: gdlat, gdlon, sm, lmh, sice
36 use gridspec_mod,
only: gridtype
37 use cmassi_mod,
only: trad_ice
38 use kinds,
only: r_kind,r_single,r_double,i_kind
39 use crtm_module,
only: crtm_atmosphere_type,crtm_surface_type,crtm_geometry_type, &
40 crtm_surface_create,o3_id,co2_id,crtm_forward,mass_mixing_ratio_units, &
41 crtm_atmosphere_create, &
42 crtm_options_type,crtm_destroy,crtm_init,specific_amount_units, &
43 success,crtm_options_destroy,crtm_options_create, crtm_options_associated
45 use crtm_rtsolution_define,
only: crtm_rtsolution_type, crtm_rtsolution_create, &
46 crtm_rtsolution_destroy, crtm_rtsolution_associated
47 use crtm_spccoeff,
only: sc
48 use crtm_atmosphere_define,
only:h2o_id,crtm_atmosphere_associated, &
49 crtm_atmosphere_destroy,volume_mixing_ratio_units,crtm_atmosphere_zero
50 use crtm_surface_define,
only: crtm_surface_destroy, crtm_surface_associated, &
52 use crtm_channelinfo_define,
only: crtm_channelinfo_type
54 use crtm_parameters,
only: limit_exp,toa_pressure,max_n_layers,max_sensor_scan_angle
55 use crtm_cloud_define,
only: water_cloud,ice_cloud,rain_cloud,snow_cloud,graupel_cloud,hail_cloud
56 use message_handler,
only: success,warning, display_message
58 use params_mod,
only: pi, rtd, p1000, capa, h1000, h1, g, rd, d608, qconv, small
59 use rqstfld_mod,
only: iget, id, lvls, iavblfld
60 use ctlblk_mod,
only: modelname, ivegsrc, novegtype, imp_physics, lm, spval, icu_physics,&
61 grib, cfld, fld_info, datapd, idat, im, jsta, jend, jm, me, ista, iend
88 INTEGER,
PARAMETER :: INVALID_LAND = 0
89 INTEGER,
PARAMETER :: COMPACTED_SOIL = 1
90 INTEGER,
PARAMETER :: TILLED_SOIL = 2
91 INTEGER,
PARAMETER :: SAND = 3
92 INTEGER,
PARAMETER :: ROCK = 4
93 INTEGER,
PARAMETER :: IRRIGATED_LOW_VEGETATION = 5
94 INTEGER,
PARAMETER :: MEADOW_GRASS = 6
95 INTEGER,
PARAMETER :: SCRUB = 7
96 INTEGER,
PARAMETER :: BROADLEAF_FOREST = 8
97 INTEGER,
PARAMETER :: PINE_FOREST = 9
98 INTEGER,
PARAMETER :: TUNDRA = 10
99 INTEGER,
PARAMETER :: GRASS_SOIL = 11
100 INTEGER,
PARAMETER :: BROADLEAF_PINE_FOREST = 12
101 INTEGER,
PARAMETER :: GRASS_SCRUB = 13
102 INTEGER,
PARAMETER :: SOIL_GRASS_SCRUB = 14
103 INTEGER,
PARAMETER :: URBAN_CONCRETE = 15
104 INTEGER,
PARAMETER :: PINE_BRUSH = 16
105 INTEGER,
PARAMETER :: BROADLEAF_BRUSH = 17
106 INTEGER,
PARAMETER :: WET_SOIL = 18
107 INTEGER,
PARAMETER :: SCRUB_SOIL = 19
108 INTEGER,
PARAMETER :: BROADLEAF70_PINE30 = 20
111 integer,
allocatable:: model_to_crtm(:)
112 integer,
parameter:: ndat=100
114 integer,
parameter:: n_absorbers = 2
116 integer,
parameter:: n_aerosols = 0
118 integer(i_kind),
parameter:: n_sensors=23
119 character(len=20),
parameter,
dimension(1:n_sensors):: sensorlist= &
143 character(len=13),
parameter,
dimension(1:n_sensors):: obslist= &
167 character(len=20),
dimension(1:n_sensors):: sensorlist_local
169 integer(i_kind) sensorindex
170 integer(i_kind) lunin,nobs,nchanl,nreal
171 integer(i_kind) error_status,itype
172 integer(i_kind) err1,err2,err3,err4
173 integer(i_kind) i,j,k,msig
174 integer(i_kind) lcbot,lctop
175 integer jdn,ichan,ixchan,igot
181 real(r_kind),
parameter:: r100=100.0_r_kind
182 real,
parameter:: ozsmall = 1.e-10
184 real(r_double),
dimension(4):: sfcpct
185 real(r_kind) snodepth,vegcover
188 real(r_kind),
dimension(ista:iend,jsta:jend):: tb1,tb2,tb3,tb4
189 real(r_kind),
allocatable :: tb(:,:,:)
190 real,
dimension(im,jm):: grid1
191 real sun_zenith,sun_azimuth, dpovg, sun_zenith_rad
194 real,
parameter:: constoz = 604229.0_r_kind
197 character(13)::obstype
199 character(20)::isis_local
201 logical hirs2,msu,goessndr,hirs3,hirs4,hirs,amsua,amsub,airs,hsb &
202 ,goes_img,abi,seviri, mhs,insat3d
203 logical avhrr,avhrr_navy,lextra,ssu
204 logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,change
205 logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img
206 logical sea,mixed,land,ice,snow,toss
207 logical micrim,microwave
208 logical post_abig16, post_abig17, post_abig18, post_abigr
209 logical fix_abig16, fix_abig17, fix_abig18
213 logical,
parameter :: debugprint = .false.
214 type(crtm_atmosphere_type),
dimension(1):: atmosphere
215 type(crtm_surface_type),
dimension(1) :: surface
216 type(crtm_geometry_type),
dimension(1) :: geometryinfo
217 type(crtm_options_type),
dimension(1) :: options
219 type(crtm_rtsolution_type),
allocatable,
dimension(:,:):: rtsolution
220 type(crtm_channelinfo_type),
allocatable,
dimension(:) :: channelinfo
222 integer ii,jj,n_clouds,n,nc
223 integer,
external :: iw3jdn
233 sensorlist_local(n) = sensorlist(n)
234 if (sensorlist(n) ==
'abi_g16')
then
235 inquire(file=
'abi_g16.SpcCoeff.bin',exist=fix_abig16)
236 if (.not.fix_abig16) sensorlist_local(n) =
'abi_gr '
238 if (sensorlist(n) ==
'abi_g17')
then
239 inquire(file=
'abi_g17.SpcCoeff.bin',exist=fix_abig17)
240 if (.not.fix_abig17) sensorlist_local(n) =
'abi_gr '
242 if (sensorlist(n) ==
'abi_g18')
then
243 inquire(file=
'abi_g18.SpcCoeff.bin',exist=fix_abig18)
244 if (.not.fix_abig18) sensorlist_local(n) =
'abi_gr '
252 allocate(model_to_crtm(novegtype) )
253 model_to_crtm=(/pine_forest, broadleaf_forest, pine_forest, &
254 broadleaf_forest,broadleaf_pine_forest, scrub, scrub_soil, &
255 broadleaf_brush,broadleaf_brush, scrub, broadleaf_brush, &
256 tilled_soil, urban_concrete,tilled_soil, urban_concrete, &
257 compacted_soil, broadleaf_brush, tundra,tundra, tundra/)
258 else if(ivegsrc==0)
then
259 allocate(model_to_crtm(novegtype) )
260 model_to_crtm=(/urban_concrete, &
261 compacted_soil, irrigated_low_vegetation, grass_soil, meadow_grass, &
262 meadow_grass, meadow_grass, scrub, grass_scrub, meadow_grass, &
263 broadleaf_forest, pine_forest, broadleaf_forest, pine_forest, &
264 broadleaf_pine_forest, compacted_soil, wet_soil, wet_soil, &
265 irrigated_low_vegetation, tundra, tundra, tundra, tundra, &
267 else if(ivegsrc==2)
then
268 allocate(model_to_crtm(0:novegtype) )
269 model_to_crtm=(/compacted_soil, &
270 broadleaf_forest, broadleaf_forest, broadleaf_pine_forest, &
271 pine_forest, pine_forest, broadleaf_brush, scrub, scrub, scrub_soil, &
272 tundra, compacted_soil, tilled_soil, compacted_soil/)
274 print*,
'novegtype=',novegtype
275 print*,
'model veg type not supported by post in calling crtm '
276 print*,
'skipping generation of simulated radiance'
284 if (iget(n) > 0) post_abig16=.true.
288 if (iget(n) > 0) post_abig17=.true.
292 if (iget(n) > 0) post_abig18=.true.
296 if (iget(n) > 0) post_abigr=.true.
300 if (iget(n) > 0) post_ahi8=.true.
304 if (iget(n) > 0) post_ssmis17=.true.
309 ifactive:
if (iget(327) > 0 .or. iget(328) > 0 .or. iget(329) > 0 &
310 .or. iget(330) > 0 .or. iget(446) > 0 .or. iget(447) > 0 &
311 .or. iget(448) > 0 .or. iget(449) > 0 .or. iget(456) > 0 &
312 .or. iget(457) > 0 .or. iget(458) > 0 .or. iget(459) > 0 &
313 .or. iget(460) > 0 .or. iget(461) > 0 .or. iget(462) > 0 &
314 .or. iget(463) > 0 .or. iget(483) > 0 .or. iget(484) > 0 &
315 .or. iget(485) > 0 .or. iget(486) > 0 .or. iget(488) > 0 &
316 .or. iget(489) > 0 .or. iget(490) > 0 .or. iget(491) > 0 &
317 .or. iget(492) > 0 .or. iget(493) > 0 .or. iget(494) > 0 &
318 .or. iget(495) > 0 .or. iget(496) > 0 .or. iget(497) > 0 &
319 .or. iget(498) > 0 .or. iget(499) > 0 .or. iget(800) > 0 &
320 .or. iget(801) > 0 .or. iget(802) > 0 .or. iget(803) > 0 &
321 .or. iget(804) > 0 .or. iget(805) > 0 .or. iget(806) > 0 &
322 .or. iget(807) > 0 .or. iget(809) > 0 &
323 .or. iget(810) > 0 .or. iget(811) > 0 .or. iget(812) > 0 &
324 .or. iget(813) > 0 .or. iget(814) > 0 .or. iget(815) > 0 &
325 .or. iget(816) > 0 .or. iget(817) > 0 .or. iget(818) > 0 &
326 .or. iget(819) > 0 .or. iget(820) > 0 .or. iget(821) > 0 &
327 .or. iget(822) > 0 .or. iget(823) > 0 .or. iget(824) > 0 &
328 .or. iget(825) > 0 .or. iget(826) > 0 .or. iget(827) > 0 &
329 .or. iget(828) > 0 .or. iget(829) > 0 .or. iget(830) > 0 &
330 .or. iget(831) > 0 .or. iget(832) > 0 .or. iget(833) > 0 &
331 .or. iget(834) > 0 .or. iget(835) > 0 .or. iget(836) > 0 &
332 .or. iget(837) > 0 .or. iget(838) > 0 .or. iget(839) > 0 &
333 .or. iget(840) > 0 .or. iget(841) > 0 .or. iget(842) > 0 &
334 .or. iget(843) > 0 .or. iget(844) > 0 .or. iget(845) > 0 &
335 .or. iget(846) > 0 .or. iget(847) > 0 .or. iget(848) > 0 &
336 .or. iget(849) > 0 .or. iget(850) > 0 .or. iget(851) > 0 &
337 .or. iget(852) > 0 .or. iget(856) > 0 .or. iget(857) > 0 &
338 .or. iget(860) > 0 .or. iget(861) > 0 &
339 .or. iget(862) > 0 .or. iget(863) > 0 .or. iget(864) > 0 &
340 .or. iget(865) > 0 .or. iget(866) > 0 .or. iget(867) > 0 &
341 .or. iget(868) > 0 .or. iget(869) > 0 .or. iget(870) > 0 &
342 .or. iget(871) > 0 .or. iget(872) > 0 .or. iget(873) > 0 &
343 .or. iget(874) > 0 .or. iget(875) > 0 .or. iget(876) > 0 &
344 .or. iget(877) > 0 .or. iget(878) > 0 .or. iget(879) > 0 &
345 .or. iget(880) > 0 .or. iget(881) > 0 .or. iget(882) > 0 &
346 .or. post_ahi8 .or. post_ssmis17 &
347 .or. post_abig16 .or. post_abig17 .or. post_abig18 &
348 .or. post_abigr )
then
350 ! specify numbers of cloud species
351 ! thompson==8, ferrier==5,95, wsm6==6, lin==2
352 if(imp_physics==99 .or. imp_physics==98)
then
354 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)
then
356 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 &
357 .or. imp_physics==28 .or. imp_physics==11)
then
361 print*,
'Warning: number of cloud species (n_clouds) being set to zero for imp_physics=',imp_physics
364 ! initialize debug print gridpoint index to middle of tile:
368 ! initialize ozone to zeros for wrf nmm and arw for now
369 if (modelname ==
'NMM' .OR. modelname ==
'NCAR' .OR. modelname ==
'RAPR' &
371 ! compute solar zenith angle for gfs, arw now computes czen in
initpost
373 jdn=iw3jdn(idat(3),idat(1),idat(2))
376 call zensun(jdn,float(idat(4)),gdlat(i,j),gdlon(i,j) &
377 ,pi,sun_zenith,sun_azimuth)
378 sun_zenith_rad=sun_zenith/rtd
379 czen(i,j)=cos(sun_zenith_rad)
382 if(jj>=jsta .and. jj<=jend.and.debugprint) &
383 print*,
'sample GFS zenith angle=',acos(czen(ii,jj))*rtd
385 ! initialize crtm. load satellite sensor array.
386 ! the
optional arguments process_id and output_process_id limit
387 ! generation of runtime informative output to mpi task
388 ! output_process_id(which here is set to be task 0)
389 if(me==0)print*,
'success in CALRAD= ',success
390 allocate( channelinfo(n_sensors))
392 error_status = crtm_init(sensorlist_local,channelinfo, &
393 process_id=0,output_process_id=0 )
394 if(me==0)print*,
'channelinfo after init= ',channelinfo(1)%sensor_id, &
395 channelinfo(2)%sensor_id
396 if (error_status /= 0_i_kind) &
397 write(6,*)
'ERROR*** crtm_init error_status=',error_status
399 ! restrict channel
list to those which are selected for simulation
400 ! in the lvls filed of wrf_cntrl.parm(does not currently apply
401 ! to all sensors / channels).
405 call select_channels_l(channelinfo(2),4,(/ 1,2,3,4 /),lvls(1:4,iget(868)),iget(868))
409 call select_channels_l(channelinfo(1),4,(/ 1,2,3,4 /),lvls(1:4,iget(872)),iget(872))
415 if (iget(n) > 0)
then
419 if (nchanl > 0 .and. nchanl <10)
then
421 if (iget(n) == 0) channelinfo(19)%Process_Channel(n-927+1)=.false.
429 if (iget(n) > 0)
then
433 if (nchanl > 0 .and. nchanl <10)
then
435 if (iget(n) == 0) channelinfo(20)%Process_Channel(n-937+1)=.false.
443 if (iget(n) > 0)
then
447 if (nchanl > 0 .and. nchanl <10)
then
449 if (iget(n) == 0) channelinfo(21)%Process_Channel(n-531+1)=.false.
453 ! goes-r for nadir output
457 if (iget(n) > 0)
then
461 if (nchanl > 0 .and. nchanl <10)
then
463 if (iget(n) == 0) channelinfo(20)%Process_Channel(n-958+1)=.false.
468 ! himawari-8 ahi infrared
472 if (iget(n) > 0)
then
476 if (nchanl > 0 .and. nchanl <10)
then
478 if (iget(n) == 0) channelinfo(22)%Process_Channel(n-969+1)=.false.
483 ! ssmis f17(37h, 37v, 91h, 91v)
487 if (iget(n) > 0)
then
491 if (nchanl > 14 .and. nchanl < 19)
then
493 if (iget(n) == 0) channelinfo(11)%Process_Channel(n-825+15)=.false.
498 ! ssmi, f13-f15(19h,19v,??h,37h,37v,85h,85v)
500 call select_channels_l(channelinfo(7),7,(/ 1,2,3,4,5,6,7 /),lvls(1:7,iget(800)),iget(800))
503 call select_channels_l(channelinfo(8),7,(/ 1,2,3,4,5,6,7 /),lvls(1:7,iget(806)),iget(806))
506 call select_channels_l(channelinfo(9),7,(/ 1,2,3,4,5,6,7 /),lvls(1:7,iget(812)),iget(812))
508 ! ssmis, f16-f20(183h,19h,19v,37h,37v,91h,91v)
510 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))
516 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))
519 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))
522 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))
526 call select_channels_l(channelinfo(15),8,(/ 1,2,3,4,5,6,7,8 /),lvls(1:8,iget(876)),iget(876))
530 call select_channels_l(channelinfo(16),4,(/ 1,2,3,4 /),lvls(1:4,iget(860)),iget(860))
534 call select_channels_l(channelinfo(17),4,(/ 1,2,3,4 /),lvls(1:4,iget(864)),iget(864))
538 call select_channels_l(channelinfo(18),4,(/ 1,2,3,4 /),lvls(1:4,iget(865)),iget(865))
541 ! loop over
data types to
process
542 sensordo:
do isat=1,n_sensors
544 if(me==0)print*,
'n_sensor,obstype,isis',isat,obslist(isat),sensorlist(isat)
546 obstype=obslist(isat)
547 isis=trim(sensorlist(isat))
550 (isis==
'imgr_g12' .and. (iget(327) > 0 .or. iget(328) > 0 &
551 .or. iget(329) > 0 .or. iget(330) > 0 .or. iget(456) > 0 &
552 .or. iget(457) > 0 .or. iget(458) > 0 .or. iget(459) > 0 )) .OR. &
553 (isis==
'imgr_g11' .and. (iget(446) > 0 .or. iget(447) > 0 &
554 .or. iget(448) > 0 .or. iget(449) > 0 .or. iget(460) > 0 &
555 .or. iget(461) > 0 .or. iget(462) > 0 .or. iget(463) > 0)) .OR. &
556 (isis==
'amsre_aqua' .and. (iget(483) > 0 .or. iget(484) > 0 &
557 .or. iget(485) > 0 .or. iget(486) > 0)) .OR. &
558 (isis==
'tmi_trmm' .and. (iget(488) > 0 .or. iget(489) > 0 &
559 .or. iget(490) > 0 .or. iget(491) > 0)) .OR. &
560 (isis==
'ssmi_f13' .and. iget(800) > 0 ) .OR. &
561 (isis==
'ssmi_f14' .and. iget(806) > 0 ) .OR. &
562 (isis==
'ssmi_f15' .and. iget(812) > 0 ) .OR. &
563 (isis==
'ssmis_f16' .and. iget(818) > 0) .OR. &
564 (isis==
'ssmis_f17' .and. post_ssmis17) .OR. &
565 (isis==
'ssmis_f18' .and. iget(832) > 0) .OR. &
566 (isis==
'ssmis_f19' .and. iget(839) > 0) .OR. &
567 (isis==
'ssmis_f20' .and. iget(846) > 0) .OR. &
568 (isis==
'imgr_mt2' .and. iget(860)>0) .OR. &
569 (isis==
'imgr_mt1r' .and. iget(864)>0) .OR. &
570 (isis==
'imgr_insat3d' .and. iget(865)>0) .OR. &
571 (isis==
'imgr_g13' .and. iget(868)>0) .OR. &
572 (isis==
'imgr_g15' .and. iget(872)>0) .OR. &
573 (isis==
'abi_g16' .and. post_abig16) .OR. &
574 (isis==
'abi_g17' .and. post_abig17) .OR. &
575 (isis==
'abi_g18' .and. post_abig18) .OR. &
576 (isis==
'abi_gr' .and. post_abigr) .OR. &
577 (isis==
'seviri_m10' .and. iget(876)>0) .OR. &
578 (isis==
'ahi_himawari8' .and. post_ahi8) )
then
579 if(me==0)print*,
'obstype, isis= ',obstype,isis
584 hirs2 = obstype ==
'hirs2'
585 hirs3 = obstype ==
'hirs3'
586 hirs4 = obstype ==
'hirs4'
587 hirs = hirs2 .or. hirs3 .or. hirs4
588 msu = obstype ==
'msu'
589 ssu = obstype ==
'ssu'
590 goessndr = obstype ==
'sndr' .or. obstype ==
'sndrd1' .or. &
591 obstype ==
'sndrd2'.or. obstype ==
'sndrd3' .or. &
593 amsua = obstype ==
'amsua'
594 amsub = obstype ==
'amsub'
595 mhs = obstype ==
'mhs'
596 airs = obstype ==
'airs'
597 hsb = obstype ==
'hsb'
598 goes_img = obstype ==
'goes_img'
599 abi = obstype ==
'abi'
600 seviri = obstype ==
'seviri'
601 insat3d = obstype ==
'imgr_insat3d'
602 avhrr = obstype ==
'avhrr'
603 avhrr_navy = obstype ==
'avhrr_navy'
604 ssmi = obstype ==
'ssmi'
605 amsre_low = obstype ==
'amsre_low'
606 amsre_mid = obstype ==
'amsre_mid'
607 amsre_hig = obstype ==
'amsre_hig'
608 amsre = amsre_low .or. amsre_mid .or. amsre_hig
609 ssmis = obstype ==
'ssmis'
610 ssmis_las = obstype ==
'ssmis_las'
611 ssmis_uas = obstype ==
'ssmis_uas'
612 ssmis_img = obstype ==
'ssmis_img'
613 ssmis_env = obstype ==
'ssmis_env'
615 ssmis=ssmis_las.or.ssmis_uas.or.ssmis_img.or.ssmis_env.or.ssmis
617 micrim=ssmi .or. ssmis .or. amsre
619 microwave=amsua .or. amsub .or. mhs .or. msu .or. hsb .or. micrim
622 sensor_search:
do j = 1, n_sensors
624 if (isis==
'abi_g16' .and. .not.fix_abig16)
then
627 if (isis==
'abi_g17' .and. .not.fix_abig17)
then
630 if (isis==
'abi_g18' .and. .not.fix_abig18)
then
633 if (channelinfo(j)%sensor_id == isis_local )
then
638 if (sensorindex == 0 )
then
639 write(6,*)
'SETUPRAD: ***WARNING*** problem with sensorindex=',isis,&
640 ' --> CAN NOT PROCESS isis=',isis,
' TERMINATE PROGRAM EXECUTION'
646 if(isis==
'ssmis_f19')channelinfo(sensorindex)%WMO_Satellite_Id=287
647 if(isis==
'ssmis_f20')channelinfo(sensorindex)%WMO_Satellite_Id=289
649 if(isis==
'abi_g16')channelinfo(sensorindex)%WMO_Satellite_Id=270
650 if(isis==
'abi_g16')channelinfo(sensorindex)%WMO_Sensor_Id=617
651 if(isis==
'abi_g17')channelinfo(sensorindex)%WMO_Satellite_Id=271
652 if(isis==
'abi_g17')channelinfo(sensorindex)%WMO_Sensor_Id=617
654 if(isis==
'abi_g18')channelinfo(sensorindex)%WMO_Satellite_Id=272
655 if(isis==
'abi_g18')channelinfo(sensorindex)%WMO_Sensor_Id=617
656 if(isis==
'abi_gr')channelinfo(sensorindex)%WMO_Satellite_Id=270
657 if(isis==
'abi_gr')channelinfo(sensorindex)%WMO_Sensor_Id=617
659 allocate(rtsolution(channelinfo(sensorindex)%n_channels,1))
660 allocate(tb(ista:iend,jsta:jend,channelinfo(sensorindex)%n_channels))
661 err1=0; err2=0; err3=0; err4=0
662 if(lm > max_n_layers)
then
663 write(6,*)
'CALRAD: lm > max_n_layers - '// &
664 'increase crtm max_n_layers ', &
669 CALL crtm_atmosphere_create(atmosphere(1),lm,n_absorbers,n_clouds &
671 CALL crtm_surface_create(surface(1),channelinfo(sensorindex)%n_channels)
672 CALL crtm_rtsolution_create(rtsolution,lm)
673 if (.NOT.(crtm_atmosphere_associated(atmosphere(1)))) &
674 write(6,*)
' ***ERROR** creating atmosphere.'
675 if (.NOT.(crtm_surface_associated(surface(1)))) &
676 write(6,*)
' ***ERROR** creating surface.'
677 if (.NOT.(any(crtm_rtsolution_associated(rtsolution)))) &
678 write(6,*)
' ***ERROR** creating rtsolution.'
680 atmosphere(1)%n_layers = lm
682 atmosphere(1)%absorber_id(1) = h2o_id
683 atmosphere(1)%absorber_id(2) = o3_id
684 atmosphere(1)%absorber_units(1) = mass_mixing_ratio_units
685 atmosphere(1)%absorber_units(2) = volume_mixing_ratio_units
687 atmosphere(1)%level_pressure(0) = toa_pressure
689 if(imp_physics==99 .or. imp_physics==98)
then
690 atmosphere(1)%cloud(1)%n_layers = lm
691 atmosphere(1)%cloud(1)%Type = water_cloud
692 atmosphere(1)%cloud(2)%n_layers = lm
693 atmosphere(1)%cloud(2)%Type = ice_cloud
694 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)
then
695 atmosphere(1)%cloud(1)%n_layers = lm
696 atmosphere(1)%cloud(1)%Type = water_cloud
697 atmosphere(1)%cloud(2)%n_layers = lm
698 atmosphere(1)%cloud(2)%Type = ice_cloud
699 atmosphere(1)%cloud(3)%n_layers = lm
700 atmosphere(1)%cloud(3)%Type = rain_cloud
701 atmosphere(1)%cloud(4)%n_layers = lm
702 atmosphere(1)%cloud(4)%Type = snow_cloud
703 atmosphere(1)%cloud(5)%n_layers = lm
704 atmosphere(1)%cloud(5)%Type = graupel_cloud
705 atmosphere(1)%cloud(6)%n_layers = lm
706 atmosphere(1)%cloud(6)%Type = hail_cloud
707 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 .or. imp_physics==28 &
708 .or. imp_physics==11)
then
709 atmosphere(1)%cloud(1)%n_layers = lm
710 atmosphere(1)%cloud(1)%Type = water_cloud
711 atmosphere(1)%cloud(2)%n_layers = lm
712 atmosphere(1)%cloud(2)%Type = ice_cloud
713 atmosphere(1)%cloud(3)%n_layers = lm
714 atmosphere(1)%cloud(3)%Type = rain_cloud
715 atmosphere(1)%cloud(4)%n_layers = lm
716 atmosphere(1)%cloud(4)%Type = snow_cloud
717 atmosphere(1)%cloud(5)%n_layers = lm
718 atmosphere(1)%cloud(5)%Type = graupel_cloud
725 surface(1)%sensordata%n_channels = channelinfo(sensorindex)%n_channels
726 surface(1)%sensordata%sensor_id = channelinfo(sensorindex)%sensor_id
727 surface(1)%sensordata%WMO_sensor_id = channelinfo(sensorindex)%WMO_sensor_id
728 surface(1)%sensordata%WMO_Satellite_id = channelinfo(sensorindex)%WMO_Satellite_id
729 surface(1)%sensordata%sensor_channel = channelinfo(sensorindex)%sensor_channel
732 nadir:
if ( (isis==
'imgr_g12' .and. (iget(327)>0 .or. &
733 iget(328)>0 .or. iget(329)>0 .or. iget(330)>0)) .or. &
734 (isis==
'imgr_g11' .and. (iget(446)>0 .or. &
735 iget(447)>0 .or. iget(448)>0 .or. iget(449)>0)) .or. &
736 (isis==
'amsre_aqua' .and. (iget(483) > 0 .or. iget(484) > 0 &
737 .or. iget(485) > 0 .or. iget(486) > 0)) .OR. &
738 (isis==
'tmi_trmm' .and. (iget(488) > 0 .or. iget(489) > 0 &
739 .or. iget(490) > 0 .or. iget(491) > 0)) .OR. &
740 (isis==
'abi_gr' .and. post_abigr) )
then
743 loopi1:
do i=ista,iend
747 if(abs(pmid(i,j,k)-spval)<=small .or. &
748 abs(t(i,j,k)-spval)<=small)
then
749 do n=1,channelinfo(sensorindex)%n_channels
758 geometryinfo(1)%sensor_zenith_angle=0.
759 geometryinfo(1)%sensor_scan_angle=0.
762 IF(geometryinfo(1)%sensor_zenith_angle <= max_sensor_scan_angle &
763 .and. geometryinfo(1)%sensor_zenith_angle >= 0.0_r_kind)
THEN
764 geometryinfo(1)%source_zenith_angle = acos(czen(i,j))*rtd
765 geometryinfo(1)%sensor_scan_angle = 0.
766 if(i==ii.and.j==jj.and.debugprint)print*,
'sample geometry ', &
767 geometryinfo(1)%sensor_zenith_angle &
768 ,geometryinfo(1)%source_zenith_angle &
771 if(modelname ==
'NCAR' .OR. modelname ==
'RAPR')
then
772 sfcpct(4)=pctsno(i,j)
773 else if(ivegsrc==1)
then
775 IF(itype == 0)itype=8
776 if(sno(i,j)<spval)
then
781 CALL snfrac (snoeqv,ivgtyp(i,j),snofrac)
783 else if(ivegsrc==2)
then
785 itype = min(max(0,ivgtyp(i,j)),13)
786 if(sno(i,j)<spval)
then
791 if(i==ii.and.j==jj.and.debugprint)print*,
'sno,itype,ivgtyp B cing snfrc = ', &
792 snoeqv,itype,ivgtyp(i,j)
793 if(sm(i,j) > 0.1)
then
796 call snfrac_gfs(snoeqv,ivgtyp(i,j),snofrac)
830 if(sm(i,j) > 0.1)
then
832 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
834 if(sfcpct(4) > 0.0_r_kind)
then
835 sfcpct(1) = 1.0_r_kind-sfcpct(4)
836 sfcpct(2) = 0.0_r_kind
837 sfcpct(3) = 0.0_r_kind
839 sfcpct(1) = 1.0_r_kind
840 sfcpct(2) = 0.0_r_kind
841 sfcpct(3) = 0.0_r_kind
844 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
846 if(sice(i,j) > 0.1)
then
847 if(sfcpct(4) > 0.0_r_kind)
then
848 sfcpct(3) = 1.0_r_kind-sfcpct(4)
849 sfcpct(1) = 0.0_r_kind
850 sfcpct(2) = 0.0_r_kind
852 sfcpct(3)= 1.0_r_kind
853 sfcpct(1) = 0.0_r_kind
854 sfcpct(2) = 0.0_r_kind
857 if(sfcpct(4) > 0.0_r_kind)
then
858 sfcpct(2)= 1.0_r_kind-sfcpct(4)
859 sfcpct(1) = 0.0_r_kind
860 sfcpct(3) = 0.0_r_kind
862 sfcpct(2)= 1.0_r_kind
863 sfcpct(1) = 0.0_r_kind
864 sfcpct(3) = 0.0_r_kind
868 if(si(i,j)/=spval)
then
875 if(ivegsrc==1 .and. itype==15 .and. sfcpct(4)<1.0_r_kind)
then
876 if(debugprint)print*,
'changing land type for veg type 15',i,j,itype,sfcpct(1:4)
884 sea = sfcpct(1) >= 0.99_r_kind
885 land = sfcpct(2) >= 0.99_r_kind
886 ice = sfcpct(3) >= 0.99_r_kind
887 snow = sfcpct(4) >= 0.99_r_kind
888 mixed = .not. sea .and. .not. ice .and. &
889 .not. land .and. .not. snow
890 if((sfcpct(1)+sfcpct(2)+sfcpct(3)+sfcpct(4)) >1._r_kind) &
891 print*,
'ERROR sfcpct ',i,j,sfcpct(1) &
892 ,sfcpct(2),sfcpct(3),sfcpct(4)
902 itype = min(max(0,ivgtyp(i,j)),novegtype)
904 itype = min(max(1,ivgtyp(i,j)),novegtype)
906 surface(1)%land_type = model_to_crtm(itype)
908 if(gridtype==
'B' .or. gridtype==
'E')
then
909 surface(1)%wind_speed = sqrt(u10h(i,j)*u10h(i,j) &
910 +v10h(i,j)*v10h(i,j))
912 surface(1)%wind_speed = sqrt(u10(i,j)*u10(i,j) &
915 surface(1)%water_coverage = sfcpct(1)
916 surface(1)%land_coverage = sfcpct(2)
917 surface(1)%ice_coverage = sfcpct(3)
918 surface(1)%snow_coverage = sfcpct(4)
920 surface(1)%land_temperature = tsfc
921 surface(1)%snow_temperature = min(tsfc,280._r_kind)
922 surface(1)%water_temperature = max(tsfc,270._r_kind)
923 surface(1)%ice_temperature = min(tsfc,280._r_kind)
924 if(smstot(i,j)/=spval)
then
925 surface(1)%soil_moisture_content = smstot(i,j)/10.
927 surface(1)%soil_moisture_content = 0.05
929 surface(1)%vegetation_fraction = vegcover
931 surface(1)%soil_temperature = 283.
933 surface(1)%snow_depth = snodepth
936 if(surface(1)%wind_speed<0. .or. surface(1)%wind_speed>200.) &
937 print*,
'bad 10 m wind'
938 if(surface(1)%water_coverage<0. .or. surface(1)%water_coverage>1.) &
939 print*,
'bad water coverage'
940 if(surface(1)%land_coverage<0. .or. surface(1)%land_coverage>1.) &
941 print*,
'bad land coverage'
942 if(surface(1)%ice_coverage<0. .or. surface(1)%ice_coverage>1.) &
943 print*,
'bad ice coverage'
944 if(surface(1)%snow_coverage<0. .or. surface(1)%snow_coverage>1.) &
945 print*,
'bad snow coverage'
946 if(surface(1)%land_temperature<0. .or. surface(1)%land_temperature>350.) &
948 if(surface(1)%soil_moisture_content<0. .or. surface(1)%soil_moisture_content>600.) &
949 print*,
'bad soil_moisture_content'
950 if(surface(1)%vegetation_fraction<0. .or. surface(1)%vegetation_fraction>1.) &
951 print*,
'bad vegetation cover'
952 if(surface(1)%snow_depth<0. .or. surface(1)%snow_depth>10000.) &
953 print*,
'bad snow_depth'
955 if(i==ii.and.j==jj.and.debugprint)print*,
'sample surface in CALRAD=', &
956 i,j,surface(1)%wind_speed,surface(1)%water_coverage, &
957 surface(1)%land_coverage,surface(1)%ice_coverage, &
958 surface(1)%snow_coverage,surface(1)%land_temperature, &
959 surface(1)%snow_temperature,surface(1)%water_temperature, &
960 surface(1)%ice_temperature,surface(1)%vegetation_fraction, &
961 surface(1)%soil_temperature,surface(1)%snow_depth, &
962 surface(1)%land_type,sm(i,j)
968 if(i==ii.and.j==jj.and.debugprint)print*,
'TOA= ',atmosphere(1)%level_pressure(0)
970 atmosphere(1)%cloud_fraction(k) = min(max(cfr(i,j,k),0.),1.)
971 atmosphere(1)%level_pressure(k) = pint(i,j,k+1)/r100
972 atmosphere(1)%pressure(k) = pmid(i,j,k)/r100
973 atmosphere(1)%temperature(k) = t(i,j,k)
974 atmosphere(1)%absorber(k,1) = max(0. &
975 ,q(i,j,k)*h1000/(h1-q(i,j,k)))
976 atmosphere(1)%absorber(k,2) = max(ozsmall,o3(i,j,k)*constoz)
980 if(atmosphere(1)%level_pressure(k)<0. .or. atmosphere(1)%level_pressure(k)>1060.) &
981 & print*,
'bad atmosphere(1)%level_pressure' &
982 & ,i,j,k,atmosphere(1)%level_pressure(k)
983 if(atmosphere(1)%pressure(k)<0. .or. &
984 & atmosphere(1)%pressure(k)>1060.) &
985 & print*,
'bad atmosphere(1)%pressure' &
986 & ,i,j,k,atmosphere(1)%pressure(k)
987 if(atmosphere(1)%temperature(k)<0. .or. &
988 & atmosphere(1)%temperature(k)>400.) &
989 & print*,
'bad atmosphere(1)%temperature'
1002 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g
1003 if(imp_physics==99 .or. imp_physics==98)
then
1004 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1005 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1009 atmosphere(1)%cloud(2)%effective_radius(k) = 50.
1010 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1012 if(atmosphere(1)%cloud(1)%water_content(k)<0. .or. &
1013 atmosphere(1)%cloud(1)%water_content(k)>1.) &
1014 print*,
'bad atmosphere cloud water'
1015 if(atmosphere(1)%cloud(2)%water_content(k)<0. .or. &
1016 atmosphere(1)%cloud(2)%water_content(k)>1.) &
1017 print*,
'bad atmosphere cloud ice'
1019 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)
then
1020 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1021 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1022 atmosphere(1)%cloud(2)%effective_radius(k) = 75.
1023 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1025 rho=pmid(i,j,k)/(rd*t(i,j,k)*(1.+d608*q(i,j,k)))
1026 atmosphere(1)%cloud(3)%effective_radius(k) = 0.
1027 if(nrain(i,j,k)>0.) &
1028 atmosphere(1)%cloud(3)%effective_radius(k) = &
1029 1.0e6*1.5*(rho*qqr(i,j,k)/(pi*rhox*nrain(i,j,k)))**(1./3.)
1030 atmosphere(1)%cloud(3)%water_content(k) = max(0.,qqr(i,j,k)*dpovg)
1031 if(f_rimef(i,j,k)<=5.0)
then
1033 if(nlice(i,j,k)>0.) &
1034 atmosphere(1)%cloud(4)%effective_radius(k) = &
1035 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1036 atmosphere(1)%cloud(4)%water_content(k) = max(0.,qqs(i,j,k)*dpovg)
1037 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1038 atmosphere(1)%cloud(5)%water_content(k) =0.
1039 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1040 atmosphere(1)%cloud(6)%water_content(k) =0.
1041 else if(f_rimef(i,j,k)<=20.0)
then
1042 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1043 atmosphere(1)%cloud(4)%water_content(k) =0.
1045 if(nlice(i,j,k)>0.) &
1046 atmosphere(1)%cloud(5)%effective_radius(k) = &
1047 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1048 atmosphere(1)%cloud(5)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1049 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1050 atmosphere(1)%cloud(6)%water_content(k) =0.
1052 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1053 atmosphere(1)%cloud(4)%water_content(k) =0.
1054 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1055 atmosphere(1)%cloud(5)%water_content(k) =0.
1057 if(nlice(i,j,k)>0.) &
1058 atmosphere(1)%cloud(6)%effective_radius(k) = &
1059 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1060 atmosphere(1)%cloud(6)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1062 if(debugprint .and. i==im/2 .and. j==jsta) &
1063 print*,
'sample precip ice radius= ',i,j,k, f_rimef(i,j,k), &
1064 atmosphere(1)%cloud(4)%effective_radius(k), atmosphere(1)%cloud(4)%water_content(k), &
1065 atmosphere(1)%cloud(5)%effective_radius(k), atmosphere(1)%cloud(5)%water_content(k), &
1066 atmosphere(1)%cloud(6)%effective_radius(k), atmosphere(1)%cloud(6)%water_content(k)
1068 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 .or. &
1069 imp_physics==28 .or. imp_physics==11)
then
1070 atmosphere(1)%cloud(1)%water_content(k)=max(0.,qqw(i,j,k)*dpovg)
1071 atmosphere(1)%cloud(2)%water_content(k)=max(0.,qqi(i,j,k)*dpovg)
1072 atmosphere(1)%cloud(3)%water_content(k)=max(0.,qqr(i,j,k)*dpovg)
1073 atmosphere(1)%cloud(4)%water_content(k)=max(0.,qqs(i,j,k)*dpovg)
1074 atmosphere(1)%cloud(5)%water_content(k)=max(0.,qqg(i,j,k)*dpovg)
1075 atmosphere(1)%cloud(1)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1076 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), &
1077 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')
1078 atmosphere(1)%cloud(2)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1079 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), &
1080 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')
1081 atmosphere(1)%cloud(3)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1082 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), &
1083 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')
1084 atmosphere(1)%cloud(4)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1085 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), &
1086 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')
1087 atmosphere(1)%cloud(5)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1088 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), &
1089 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')
1094 if (abs(atmosphere(1)%pressure(k)-atmosphere(1)%pressure(k+1)) &
1096 atmosphere(1)%pressure(k)=atmosphere(1)%pressure(k)+0.005
1102 if(icu_physics==2)
then
1103 lcbot=nint(hbot(i,j))
1104 lctop=nint(htop(i,j))
1105 if (lcbot-lctop > 1)
then
1109 q_conv = cnvcfr(i,j)*qconv
1111 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g
1112 if (t(i,j,k) < trad_ice)
then
1113 atmosphere(1)%cloud(2)%water_content(k) = &
1114 atmosphere(1)%cloud(2)%water_content(k) + dpovg*q_conv
1116 atmosphere(1)%cloud(1)%water_content(k) = &
1117 atmosphere(1)%cloud(1)%water_content(k) + dpovg*q_conv
1125 error_status = crtm_forward(atmosphere,surface, &
1126 geometryinfo,channelinfo(sensorindex:sensorindex), &
1128 if (error_status /=0)
then
1129 print*,
'***ERROR*** during crtm_forward call ', error_status
1130 do n=1,channelinfo(sensorindex)%n_channels
1137 do n=1,channelinfo(sensorindex)%n_channels
1138 tb(i,j,n)=rtsolution(n,1)%brightness_temperature
1144 if(i==ii.and.j==jj)
then
1145 do n=1,channelinfo(sensorindex)%n_channels
1146 3301
format(
'Sample rtsolution(',i0,
',',i0,
') in CALRAD = ',f0.3)
1147 print 3301,n,1,rtsolution(n,1)%brightness_temperature
1149 do n=1,channelinfo(sensorindex)%n_channels
1150 3302
format(
'Sample tb(',i0,
',',i0,
',',i0,
') in CALRAD = ',f0.3)
1151 print 3302,ii,jj,n,tb(ii,jj,n)
1161 do n=1,channelinfo(sensorindex)%n_channels
1176 if (isis==
'amsre_aqua')
then
1179 igot=iget(482+ixchan)
1183 grid1(i,j)=tb(i,j,ichan)
1186 if (grib==
"grib2")
then
1188 fld_info(cfld)%ifld=iavblfld(igot)
1189 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1194 if (isis==
'tmi_trmm')
then
1197 igot=iget(487+ixchan)
1201 grid1(i,j) = tb(i,j,ichan)
1204 if (grib==
"grib2")
then
1206 fld_info(cfld)%ifld=iavblfld(igot)
1207 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1213 if (isis==
'imgr_g11')
then
1220 grid1(i,j) = tb(i,j,ichan)
1223 if (grib==
"grib2")
then
1225 fld_info(cfld)%ifld=iavblfld(igot)
1226 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1232 if (isis==
'imgr_g12')
then
1235 igot=iget(326+ixchan)
1239 grid1(i,j)=tb(i,j,ichan)
1242 if (grib==
"grib2")
then
1244 fld_info(cfld)%ifld=iavblfld(igot)
1245 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1250 if (isis==
'abi_gr')
then
1254 igot=iget(957+ixchan)
1258 grid1(i,j)=tb(i,j,ichan)
1261 if(grib==
"grib2" )
then
1263 fld_info(cfld)%ifld=iavblfld(igot)
1264 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1275 nonnadir:
if((isis==
'ssmi_f13' .and. iget(800) > 0 ) .OR. &
1276 (isis==
'ssmi_f14' .and. iget(806) > 0 ) .OR. &
1277 (isis==
'ssmi_f15' .and. iget(812) > 0 ) .OR. &
1278 (isis==
'ssmis_f16' .and. iget(818) > 0) .OR. &
1279 (isis==
'ssmis_f17' .and. post_ssmis17) .OR. &
1280 (isis==
'ssmis_f18' .and. iget(832) > 0) .OR. &
1281 (isis==
'ssmis_f19' .and. iget(839) > 0) .OR. &
1282 (isis==
'ssmis_f20' .and. iget(846) > 0) .OR. &
1283 (isis==
'imgr_mt2' .and. iget(860)>0) .OR. &
1284 (isis==
'imgr_mt1r' .and. iget(864)>0) .OR. &
1285 (isis==
'imgr_insat3d' .and. iget(865)>0) .OR. &
1286 (isis==
'imgr_g13' .and. iget(868)>0) .OR. &
1287 (isis==
'imgr_g15' .and. iget(872)>0) .OR. &
1288 (isis==
'abi_g16' .and. post_abig16) .OR. &
1289 (isis==
'abi_g17' .and. post_abig17) .OR. &
1290 (isis==
'abi_g18' .and. post_abig18) .OR. &
1291 (isis==
'seviri_m10' .and. iget(876)>0) .OR. &
1292 (isis==
'ahi_himawari8' .and. post_ahi8) .OR. &
1293 (isis==
'imgr_g12' .and. (iget(456)>0 .or. &
1294 iget(457)>0 .or. iget(458)>0 .or. iget(459)>0)) .or. &
1295 (isis==
'imgr_g11' .and. (iget(460)>0 .or. &
1296 iget(461)>0 .or. iget(462)>0 .or. iget(463)>0)))
then
1299 loopi2:
do i=ista,iend
1303 if(abs(pmid(i,j,k)-spval)<=small .or. &
1304 abs(t(i,j,k)-spval)<=small)
then
1305 do n=1,channelinfo(sensorindex)%n_channels
1315 if(isis==
'imgr_g12')
then
1318 else if(isis==
'seviri_m10')
then
1321 else if(isis==
'imgr_g13')
then
1324 else if(isis==
'imgr_g15')
then
1327 else if(isis==
'abi_g16')
then
1330 else if(isis==
'abi_g17')
then
1333 else if(isis==
'abi_g18')
then
1336 else if(isis==
'imgr_g11')
then
1339 else if(isis==
'imgr_mt2')
then
1342 else if(isis==
'imgr_mt1r')
then
1345 else if(isis==
'imgr_insat3d')
then
1348 else if(isis==
'ahi_himawari8')
then
1354 if(isis==
'ssmis_f16'.or.isis==
'ssmis_f17'.or.isis==
'ssmis_f18'.or. &
1355 isis==
'ssmis_f19'.or.isis==
'ssmis_f20'.or.isis==
'ssmi_f13'.or. &
1356 isis==
'ssmi_f14'.or.isis==
'ssmi_f15')
then
1360 call geo_zenith_angle(i,j,gdlat(i,j),gdlon(i,j) &
1361 ,sublat,sublon,sat_zenith)
1364 geometryinfo(1)%sensor_zenith_angle=sat_zenith
1365 geometryinfo(1)%sensor_scan_angle=sat_zenith
1367 if(i==ii .and. j==jj.and.debugprint)
then
1368 print *,
'zenith info: zenith=',sat_zenith,
' scan=',sat_zenith, &
1369 ' MAX_SENSOR_SCAN_ANGLE=',max_sensor_scan_angle
1373 IF(geometryinfo(1)%sensor_zenith_angle <= max_sensor_scan_angle &
1374 .and. geometryinfo(1)%sensor_zenith_angle >= 0.0_r_kind)
THEN
1375 geometryinfo(1)%source_zenith_angle = acos(czen(i,j))*rtd
1376 geometryinfo(1)%sensor_scan_angle = 0.
1377 if(i==ii.and.j==jj.and.debugprint)print*,
'sample geometry ', &
1378 geometryinfo(1)%sensor_zenith_angle &
1379 ,geometryinfo(1)%source_zenith_angle &
1383 if(modelname ==
'NCAR' .OR. modelname ==
'RAPR')
then
1384 sfcpct(4)=pctsno(i,j)
1385 else if(ivegsrc==1)
then
1387 IF(itype == 0)itype=8
1388 if(sno(i,j)<spval)
then
1393 CALL snfrac (snoeqv,ivgtyp(i,j),snofrac)
1395 else if(ivegsrc==2)
then
1397 itype = min(max(0,ivgtyp(i,j)),13)
1398 if(sno(i,j)<spval)
then
1403 if(i==ii.and.j==jj)print*,
'sno,itype,ivgtyp B cing snfrc = ', &
1404 snoeqv,itype,ivgtyp(i,j)
1405 if(sm(i,j) > 0.1)
then
1408 call snfrac_gfs(snoeqv,ivgtyp(i,j),snofrac)
1415 if(sm(i,j) > 0.1)
then
1417 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
1419 if(sfcpct(4) > 0.0_r_kind)
then
1420 sfcpct(1) = 1.0_r_kind-sfcpct(4)
1421 sfcpct(2) = 0.0_r_kind
1422 sfcpct(3) = 0.0_r_kind
1424 sfcpct(1) = 1.0_r_kind
1425 sfcpct(2) = 0.0_r_kind
1426 sfcpct(3) = 0.0_r_kind
1429 tsfc = ths(i,j)*(pint(i,j,nint(lmh(i,j))+1)/p1000)**capa
1430 vegcover=vegfrc(i,j)
1431 if(sice(i,j) > 0.1)
then
1432 if(sfcpct(4) > 0.0_r_kind)
then
1433 sfcpct(3) = 1.0_r_kind-sfcpct(4)
1434 sfcpct(1) = 0.0_r_kind
1435 sfcpct(2) = 0.0_r_kind
1437 sfcpct(3)= 1.0_r_kind
1438 sfcpct(1) = 0.0_r_kind
1439 sfcpct(2) = 0.0_r_kind
1442 if(sfcpct(4) > 0.0_r_kind)
then
1443 sfcpct(2)= 1.0_r_kind-sfcpct(4)
1444 sfcpct(1) = 0.0_r_kind
1445 sfcpct(3) = 0.0_r_kind
1447 sfcpct(2)= 1.0_r_kind
1448 sfcpct(1) = 0.0_r_kind
1449 sfcpct(3) = 0.0_r_kind
1453 if(si(i,j)/=spval)
then
1462 if(novegtype==20 .and. itype==15 .and.sfcpct(4)<1.0_r_kind)
then
1463 if(debugprint)print*,
'changing land type for veg type 15',i,j,itype,sfcpct(1:4)
1464 sfcpct(1)=0.0_r_kind
1465 sfcpct(2)=0.0_r_kind
1466 sfcpct(3)=0.0_r_kind
1467 sfcpct(4)=1.0_r_kind
1472 sea = sfcpct(1) >= 0.99_r_kind
1473 land = sfcpct(2) >= 0.99_r_kind
1474 ice = sfcpct(3) >= 0.99_r_kind
1475 snow = sfcpct(4) >= 0.99_r_kind
1476 mixed = .not. sea .and. .not. ice .and. &
1477 .not. land .and. .not. snow
1478 if((sfcpct(1)+sfcpct(2)+sfcpct(3)+sfcpct(4)) >1._r_kind) &
1479 print*,
'ERROR sfcpct ',i,j,sfcpct(1) &
1480 ,sfcpct(2),sfcpct(3),sfcpct(4)
1490 itype = min(max(0,ivgtyp(i,j)),novegtype)
1492 itype = min(max(1,ivgtyp(i,j)),novegtype)
1494 surface(1)%land_type = model_to_crtm(itype)
1496 if(gridtype==
'B' .or. gridtype==
'E')
then
1497 surface(1)%wind_speed = sqrt(u10h(i,j)*u10h(i,j) &
1498 +v10h(i,j)*v10h(i,j))
1500 surface(1)%wind_speed = sqrt(u10(i,j)*u10(i,j) &
1504 surface(1)%water_coverage = sfcpct(1)
1505 surface(1)%land_coverage = sfcpct(2)
1506 surface(1)%ice_coverage = sfcpct(3)
1507 surface(1)%snow_coverage = sfcpct(4)
1508 surface(1)%land_temperature = tsfc
1509 surface(1)%snow_temperature = min(tsfc,280._r_kind)
1510 surface(1)%water_temperature = max(tsfc,270._r_kind)
1511 surface(1)%ice_temperature = min(tsfc,280._r_kind)
1512 if(smstot(i,j)/=spval)
then
1513 surface(1)%soil_moisture_content = smstot(i,j)/10.
1515 surface(1)%soil_moisture_content = 0.05
1517 surface(1)%vegetation_fraction = vegcover
1519 surface(1)%soil_temperature = 283.
1521 surface(1)%snow_depth = snodepth
1524 if(surface(1)%wind_speed<0. .or. surface(1)%wind_speed>200.) &
1525 print*,
'bad 10 m wind'
1526 if(surface(1)%water_coverage<0. .or. surface(1)%water_coverage>1.) &
1527 print*,
'bad water coverage'
1528 if(surface(1)%land_coverage<0. .or. surface(1)%land_coverage>1.) &
1529 print*,
'bad land coverage'
1530 if(surface(1)%ice_coverage<0. .or. surface(1)%ice_coverage>1.) &
1531 print*,
'bad ice coverage'
1532 if(surface(1)%snow_coverage<0. .or. surface(1)%snow_coverage>1.) &
1533 print*,
'bad snow coverage'
1534 if(surface(1)%land_temperature<0. .or. surface(1)%land_temperature>350.) &
1536 if(surface(1)%soil_moisture_content<0. .or. surface(1)%soil_moisture_content>600.) &
1537 print*,
'bad soil_moisture_content'
1538 if(surface(1)%vegetation_fraction<0. .or. surface(1)%vegetation_fraction>1.) &
1539 print*,
'bad vegetation cover'
1540 if(surface(1)%snow_depth<0. .or. surface(1)%snow_depth>10000.) &
1541 print*,
'bad snow_depth'
1544 if(i==ii.and.j==jj.and.debugprint)print*,
'sample surface in CALRAD=', &
1545 i,j,surface(1)%wind_speed,surface(1)%water_coverage, &
1546 surface(1)%land_coverage,surface(1)%ice_coverage, &
1547 surface(1)%snow_coverage,surface(1)%land_temperature, &
1548 surface(1)%snow_temperature,surface(1)%water_temperature, &
1549 surface(1)%ice_temperature,surface(1)%vegetation_fraction, &
1550 surface(1)%soil_temperature,surface(1)%snow_depth, &
1551 surface(1)%land_type,sm(i,j)
1557 if(i==ii.and.j==jj.and.debugprint)print*,
'TOA= ',atmosphere(1)%level_pressure(0)
1559 atmosphere(1)%cloud_fraction(k) = min(max(cfr(i,j,k),0.),1.)
1560 atmosphere(1)%level_pressure(k) = pint(i,j,k+1)/r100
1561 atmosphere(1)%pressure(k) = pmid(i,j,k)/r100
1562 atmosphere(1)%temperature(k) = t(i,j,k)
1563 atmosphere(1)%absorber(k,1) = max(0. &
1564 ,q(i,j,k)*h1000/(h1-q(i,j,k)))
1565 atmosphere(1)%absorber(k,2) = max(ozsmall,o3(i,j,k)*constoz)
1569 if(atmosphere(1)%level_pressure(k)<0. .or. atmosphere(1)%level_pressure(k)>1060.) &
1570 print*,
'bad atmosphere(1)%level_pressure' &
1571 ,i,j,k,atmosphere(1)%level_pressure(k)
1572 if(atmosphere(1)%pressure(k)<0. .or. &
1573 atmosphere(1)%pressure(k)>1060.) &
1574 print*,
'bad atmosphere(1)%pressure' &
1575 ,i,j,k,atmosphere(1)%pressure(k)
1576 if(atmosphere(1)%temperature(k)<0. .or. &
1577 atmosphere(1)%temperature(k)>400.) &
1578 print*,
'bad atmosphere(1)%temperature'
1591 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g
1592 if(imp_physics==99 .or. imp_physics==98)
then
1593 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1594 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1598 atmosphere(1)%cloud(2)%effective_radius(k) = 50.
1599 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1601 if(atmosphere(1)%cloud(1)%water_content(k)<0. .or. &
1602 atmosphere(1)%cloud(1)%water_content(k)>1.) &
1603 print*,
'bad atmosphere cloud water'
1604 if(atmosphere(1)%cloud(2)%water_content(k)<0. .or. &
1605 atmosphere(1)%cloud(2)%water_content(k)>1.) &
1606 print*,
'bad atmosphere cloud ice'
1608 else if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)
then
1609 atmosphere(1)%cloud(1)%effective_radius(k) = 10.
1610 atmosphere(1)%cloud(1)%water_content(k) = max(0.,qqw(i,j,k)*dpovg)
1611 atmosphere(1)%cloud(2)%effective_radius(k) = 75.
1612 atmosphere(1)%cloud(2)%water_content(k) = max(0.,qqi(i,j,k)*dpovg)
1614 rho=pmid(i,j,k)/(rd*t(i,j,k)*(1.+d608*q(i,j,k)))
1615 atmosphere(1)%cloud(3)%effective_radius(k) = 0.
1616 if(nrain(i,j,k)>0.) &
1617 atmosphere(1)%cloud(3)%effective_radius(k) = &
1618 1.0e6*1.5*(rho*qqr(i,j,k)/(pi*rhox*nrain(i,j,k)))**(1./3.)
1619 atmosphere(1)%cloud(3)%water_content(k) = max(0.,qqr(i,j,k)*dpovg)
1620 if(f_rimef(i,j,k)<=5.0)
then
1622 if(nlice(i,j,k)>0.) &
1623 atmosphere(1)%cloud(4)%effective_radius(k) = &
1624 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1625 atmosphere(1)%cloud(4)%water_content(k) = max(0.,qqs(i,j,k)*dpovg)
1626 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1627 atmosphere(1)%cloud(5)%water_content(k) =0.
1628 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1629 atmosphere(1)%cloud(6)%water_content(k) =0.
1630 else if(f_rimef(i,j,k)<=20.0)
then
1631 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1632 atmosphere(1)%cloud(4)%water_content(k) =0.
1634 if(nlice(i,j,k)>0.) &
1635 atmosphere(1)%cloud(5)%effective_radius(k) = &
1636 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1637 atmosphere(1)%cloud(5)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1638 atmosphere(1)%cloud(6)%effective_radius(k) = 0.
1639 atmosphere(1)%cloud(6)%water_content(k) =0.
1641 atmosphere(1)%cloud(4)%effective_radius(k) = 0.
1642 atmosphere(1)%cloud(4)%water_content(k) =0.
1643 atmosphere(1)%cloud(5)%effective_radius(k) = 0.
1644 atmosphere(1)%cloud(5)%water_content(k) =0.
1646 if(nlice(i,j,k)>0.) &
1647 atmosphere(1)%cloud(6)%effective_radius(k) = &
1648 1.0e6*1.5*(rho*qqs(i,j,k)/(pi*rhox*nlice(i,j,k)))**(1./3.)
1649 atmosphere(1)%cloud(6)%water_content(k) =max(0.,qqs(i,j,k)*dpovg)
1651 if(debugprint .and. i==im/2 .and. j==jsta) &
1652 print*,
'sample precip ice radius= ',i,j,k, f_rimef(i,j,k), &
1653 atmosphere(1)%cloud(4)%effective_radius(k), atmosphere(1)%cloud(4)%water_content(k), &
1654 atmosphere(1)%cloud(5)%effective_radius(k), atmosphere(1)%cloud(5)%water_content(k), &
1655 atmosphere(1)%cloud(6)%effective_radius(k), atmosphere(1)%cloud(6)%water_content(k)
1657 else if(imp_physics==8 .or. imp_physics==6 .or. imp_physics==2 .or. &
1658 imp_physics==28 .or. imp_physics==11)
then
1659 atmosphere(1)%cloud(1)%water_content(k)=max(0.,qqw(i,j,k)*dpovg)
1660 atmosphere(1)%cloud(2)%water_content(k)=max(0.,qqi(i,j,k)*dpovg)
1661 atmosphere(1)%cloud(3)%water_content(k)=max(0.,qqr(i,j,k)*dpovg)
1662 atmosphere(1)%cloud(4)%water_content(k)=max(0.,qqs(i,j,k)*dpovg)
1663 atmosphere(1)%cloud(5)%water_content(k)=max(0.,qqg(i,j,k)*dpovg)
1665 if(effrl(i,j,k)/=spval)
then
1666 atmosphere(1)%cloud(1)%effective_radius(k)=min(max(effrl(i,j,k),2.5),50.)
1668 atmosphere(1)%cloud(1)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1669 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), &
1670 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')
1672 if(effri(i,j,k)/=spval)
then
1673 atmosphere(1)%cloud(2)%effective_radius(k)=min(max(effri(i,j,k),2.5),125.)
1675 atmosphere(1)%cloud(2)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1676 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), &
1677 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')
1679 if(effrs(i,j,k)/=spval)
then
1680 atmosphere(1)%cloud(4)%effective_radius(k)=min(max(effrs(i,j,k),2.5),1000.)
1682 atmosphere(1)%cloud(4)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1683 q(i,j,k),qqw(i,j,k),qqi(i,j,k),qqr(i,j,k),f_rimef(i,j,k),nlice(i,j,k), &
1684 nrain(i,j,k),qqs(i,j,k),qqg(i,j,k),qqnr(i,j,k),qqni(i,j,k),qqnw(i,j,k),imp_physics,
'S')
1686 atmosphere(1)%cloud(3)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1687 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), &
1688 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')
1689 atmosphere(1)%cloud(5)%effective_radius(k)=effr(pmid(i,j,k),t(i,j,k), &
1690 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), &
1691 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')
1696 if (abs(atmosphere(1)%pressure(k)-atmosphere(1)%pressure(k+1)) &
1698 atmosphere(1)%pressure(k)=atmosphere(1)%pressure(k)+0.005
1703 if(icu_physics==2)
then
1704 lcbot=nint(hbot(i,j))
1705 lctop=nint(htop(i,j))
1706 if (lcbot-lctop > 1)
then
1710 q_conv = cnvcfr(i,j)*qconv
1712 dpovg=(pint(i,j,k+1)-pint(i,j,k))/g
1713 if (t(i,j,k) < trad_ice)
then
1714 atmosphere(1)%cloud(2)%water_content(k) = &
1715 atmosphere(1)%cloud(2)%water_content(k) + dpovg*q_conv
1717 atmosphere(1)%cloud(1)%water_content(k) = &
1718 atmosphere(1)%cloud(1)%water_content(k) + dpovg*q_conv
1726 error_status = crtm_forward(atmosphere,surface, &
1727 geometryinfo,channelinfo(sensorindex:sensorindex), &
1729 if (error_status /=0)
then
1730 print*,
'***ERROR*** during crtm_forward call ', &
1732 do n=1,channelinfo(sensorindex)%n_channels
1736 do n=1,channelinfo(sensorindex)%n_channels
1737 tb(i,j,n)=rtsolution(n,1)%brightness_temperature
1739 if(i==ii.and.j==jj.and.debugprint)
then
1740 do n=1,channelinfo(sensorindex)%n_channels
1741 3303
format(
'Sample rtsolution(',i0,
',',i0,
') in CALRAD = ',f0.3)
1744 do n=1,channelinfo(sensorindex)%n_channels
1745 3304
format(
'Sample tb(',i0,
',',i0,
',',i0,
') in CALRAD = ',f0.3)
1746 print 3304,i,j,n,tb(i,j,n)
1756 do n=1,channelinfo(sensorindex)%n_channels
1767 if (isis==
'ssmi_f13')
then
1773 if(lvls(ixchan,igot)==1)
then
1777 grid1(i,j)=tb(i,j,nc)
1780 if (grib==
"grib2")
then
1782 fld_info(cfld)%ifld=iavblfld(igot)
1783 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1789 if (isis==
'ssmi_f14')
then
1795 if(lvls(ixchan,igot)==1)
then
1799 grid1(i,j)=tb(i,j,nc)
1802 if (grib==
"grib2")
then
1804 fld_info(cfld)%ifld=iavblfld(igot)
1805 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1811 if (isis==
'ssmi_f15')
then
1817 if(lvls(ixchan,igot)==1)
then
1821 grid1(i,j)=tb(i,j,nc)
1824 if (grib==
"grib2")
then
1826 fld_info(cfld)%ifld=iavblfld(igot)
1827 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1833 if (isis==
'ssmis_f16')
then
1839 print*,
'ixchan,lvls=',ixchan,lvls(ixchan,igot)
1840 if(lvls(ixchan,igot)==1)
then
1844 grid1(i,j)=tb(i,j,nc)
1847 if (grib==
"grib2")
then
1849 fld_info(cfld)%ifld=iavblfld(igot)
1850 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1856 if(isis==
'ssmis_f17')
then
1859 igot=iget(824+ixchan)
1863 grid1(i,j)=tb(i,j,ichan)
1866 if(grib==
"grib2" )
then
1868 fld_info(cfld)%ifld=iavblfld(igot)
1869 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1874 if (isis==
'ssmis_f18')
then
1880 if(lvls(ixchan,igot)==1)
then
1884 grid1(i,j)=tb(i,j,nc)
1887 if (grib==
"grib2")
then
1889 fld_info(cfld)%ifld=iavblfld(igot)
1890 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1896 if (isis==
'ssmis_f19')
then
1902 if(lvls(ixchan,igot)==1)
then
1906 grid1(i,j)=tb(i,j,nc)
1909 if (grib==
"grib2")
then
1911 fld_info(cfld)%ifld=iavblfld(igot)
1912 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1918 if (isis==
'ssmis_f20')
then
1924 if(lvls(ixchan,igot)==1)
then
1928 grid1(i,j)=tb(i,j,nc)
1931 if (grib==
"grib2")
then
1933 fld_info(cfld)%ifld=iavblfld(igot)
1934 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1940 if(isis==
'imgr_mt2')
then
1944 if(lvls(ichan,igot)==1)
then
1948 grid1(i,j)=tb(i,j,nc)
1951 if(grib==
"grib2")
then
1953 fld_info(cfld)%ifld=iavblfld(igot)
1954 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1959 if(isis==
'imgr_mt1r')
then
1963 if(lvls(ichan,igot)==1)
then
1967 grid1(i,j)=tb(i,j,nc)
1970 if(grib==
"grib2" )
then
1972 fld_info(cfld)%ifld=iavblfld(igot)
1973 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1978 if_insat3d:
if(isis==
'imgr_insat3d')
then
1982 if(lvls(ichan,igot)==1)
then
1986 grid1(i,j)=tb(i,j,nc)
1989 if(grib==
"grib2" )
then
1991 fld_info(cfld)%ifld=iavblfld(igot)
1992 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1997 if (isis==
'imgr_g11')
then
2000 igot=iget(459+ixchan)
2004 grid1(i,j)=tb(i,j,ichan)
2007 if(grib==
"grib2" )
then
2009 fld_info(cfld)%ifld=iavblfld(igot)
2010 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2015 if (isis==
'imgr_g12')
then
2018 igot=iget(455+ixchan)
2022 grid1(i,j)=tb(i,j,ichan)
2025 if(grib==
"grib2" )
then
2027 fld_info(cfld)%ifld=iavblfld(igot)
2028 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2033 if (isis==
'seviri_m10')
then
2039 if(lvls(ixchan,igot)==1)
then
2043 grid1(i,j)=tb(i,j,nc)
2046 if (grib==
"grib2")
then
2048 fld_info(cfld)%ifld=iavblfld(igot)
2049 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2055 if (isis==
'imgr_g13')
then
2061 if(lvls(ixchan,igot)==1)
then
2065 grid1(i,j)=tb(i,j,nc)
2068 if (grib==
"grib2")
then
2070 fld_info(cfld)%ifld=iavblfld(igot)
2071 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2077 if (isis==
'imgr_g15')
then
2083 if(lvls(ixchan,igot)==1)
then
2087 grid1(i,j)=tb(i,j,nc)
2090 if (grib==
"grib2")
then
2092 fld_info(cfld)%ifld=iavblfld(igot)
2093 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2099 if (isis==
'abi_g16')
then
2103 igot=iget(926+ixchan)
2107 grid1(i,j)=tb(i,j,ichan)
2110 if(grib==
"grib2" )
then
2112 fld_info(cfld)%ifld=iavblfld(igot)
2113 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2118 if (isis==
'abi_g17')
then
2122 igot=iget(936+ixchan)
2126 grid1(i,j)=tb(i,j,ichan)
2129 if(grib==
"grib2" )
then
2131 fld_info(cfld)%ifld=iavblfld(igot)
2132 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2138 if (isis==
'abi_g18')
then
2141 igot=iget(530+ixchan)
2146 grid1(i,j)=tb(i,j,ichan)
2153 if(grib==
"grib2" )
then
2155 fld_info(cfld)%ifld=iavblfld(igot)
2156 datapd(1:im,1:jend-jsta+1,cfld)=grid1(1:im,jsta:jend)
2161 if(isis==
'ahi_himawari8')
then
2163 igot=iget(968+ichan)
2167 grid1(i,j)=tb(i,j,ichan)
2170 if(grib==
"grib2" )
then
2172 fld_info(cfld)%ifld=iavblfld(igot)
2173 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2182 CALL crtm_atmosphere_destroy(atmosphere(1))
2183 CALL crtm_surface_destroy(surface(1))
2184 CALL crtm_rtsolution_destroy(rtsolution)
2185 if (crtm_atmosphere_associated(atmosphere(1))) &
2186 write(6,*)
' ***ERROR** destroying atmosphere.'
2187 if (crtm_surface_associated(surface(1))) &
2188 write(6,*)
' ***ERROR** destroying surface.'
2189 if (any(crtm_rtsolution_associated(rtsolution))) &
2190 write(6,*)
' ***ERROR** destroying rtsolution.'
2192 deallocate (rtsolution)
2197 error_status = crtm_destroy(channelinfo)
2198 if (error_status /= success) &
2199 print*,
'ERROR*** crtm_destroy error_status=',error_status
2200 deallocate(channelinfo)
2201 if (
allocated(model_to_crtm))
deallocate(model_to_crtm)