1module multiIRbgck_mod
2 ! MODULE multiIRbgck_mod (prefix='irbg' category='1. High-level functionality')
3 !
4 !:Purpose: Variables for multispectral infrared background check and quality
5 ! control.
6 !
7 use rttovInterfaces_mod
8 use tovsNL_mod
9 use rttov_const, only : inst_id_iasi
10 use rttov_types, only : &
11 rttov_coefs ,&
12 rttov_profile ,&
13 rttov_radiance ,&
14 rttov_transmission ,&
15 rttov_chanprof ,&
16 rttov_emissivity
17 use utilities_mod
18 use obsSpaceData_mod
19 use midasMpi_mod
20 use columnData_mod
21 use MathPhysConstants_mod
22 implicit none
23 save
24 private
25 ! Public functions (methods)
26 public :: irbg_setup, irbg_bgCheckIR
27
28 integer, parameter :: nClassAVHRR = 7
29 integer, parameter :: NIR = 3, NVIS = 3
30 integer, parameter :: nChanAVHRR = NIR + NVIS
31 integer, parameter :: nmaxinst = 10
32
33 ! Cloud top units : (1) mb, (2) meters
34 ! (subroutines cloud_height (iopt1) and cloud_top (iopt2))
35
36 integer, parameter :: iopt1 = 2 ! verify subr input if iopt1 changes
37 integer, parameter :: iopt2 = 1
38
39 ! Cloud top based on which background profile matching (subroutine cloud_top)
40 ! (0) brightness temperature, (1) radiance, (2) both
41
42 integer, parameter :: ihgt = 2
43
44 ! Highest flag in post files (value of N in 2^N)
45
46 integer, parameter :: bitflag = 29
47
48 real(8),parameter :: albedoThresholdQC= 0.25d0 ! max albedo allowed over the water
49
50 real(8),parameter :: seuilalb_static(NIR,0:2)= reshape( (/ 70.0,67.0,50.0, &
51 40.0,37.0,37.0, &
52 70.0,57.0,40. /),(/3,3/) )
53 real(8),parameter :: seuilalb_homog(NIR,0:2)= reshape( (/ 15.0,18.0,13.0, &
54 9.0,10.0,10.0, &
55 18.0,16.0,10.0 /),(/3,3/) )
56 real(8) :: seuilbt_homog(NVIS+1:nChanAVHRR,0:2,1:2)= reshape( (/5.d0, 4.d0, 4.d0, 4.d0, 3.d0, 3.d0, &
57 5.d0, 4.d0, 4.d0, 5.d0, 5.d0, 5.d0, &
58 4.d0, 3.d0, 3.d0, 5.d0, 5.d0, 5.d0/), (/3,3,2/) )
59
60 ! Number of channels to use for cloud top height detection
61 ! with the "background profile matching" method (subroutine cloud_top)
62 integer, parameter :: nch_he = 4
63
64 ! Number of channels pairs to use for cloud top height detection
65 ! with the CO2-slicing method. (subroutine co2_slicing)
66 integer, parameter :: nco2 = 13
67
68 ! Namelist variables:
69 integer :: ninst ! MUST NOT BE INCLUDED IN NAMELIST!
70 character(len=7) :: inst(nmaxinst) ! List of instrument names
71 integer :: iwindow(nmaxinst) ! Ref window channel for clear/cloudy profile detection
72 integer :: iwindow_alt(nmaxinst) ! Alternate window channel for clear/cloudy profile detection
73 integer :: ilist1(nmaxinst,nch_he) ! Chan numbers for cloud top height detection: background profile matching
74 integer :: ilist2(nmaxinst,nco2) ! Chan numbers for cloud top height detection: CO2-slicing
75 integer :: ilist2_pair(nmaxinst,nco2) ! Chan number pairs for cloud top height detection: CO2-slicing
76 real(8) :: dtw ! Max delta allowed btwn guess and true skin temp over water
77 real(8) :: dtl ! Max delta allowed btwn guess and true skin temp over land
78 real(8) :: pco2min ! Min RTTOV level for lev_start variable entering CO2 slicing in mb
79 real(8) :: pco2max ! Max RTTOV level for lev_start variable entering CO2 slicing in mb
80 real(8) :: night_ang ! Min solar zenith angle for night (between 90 and 180)
81 real(8) :: crisCloudFractionThreshold ! threshold for CrIS cloud detection from VIIRS cloud mask
82
83 type(rttov_coefs) :: coefs_avhrr
84
85 type avhrr_bgck_iasi
86 real(8) :: radmoy(nClassAVHRR,nChanAVHRR)
87 real(8) :: radstd(nClassAVHRR,nChanAVHRR)
88 real(8) :: cfrac(nClassAVHRR)
89 real(8) :: tbmoy(nClassAVHRR,NVIS+1:NVIS+NIR)
90 real(8) :: tbstd(nClassAVHRR,NVIS+1:NVIS+NIR)
91 real(8) :: albedmoy(nClassAVHRR,1:NVIS)
92 real(8) :: albedstd(nClassAVHRR,1:NVIS)
93 real(8) :: tbstd_pixelIASI(NVIS+1:NVIS+NIR)
94 real(8) :: albstd_pixeliasi(1:NVIS)
95 real(8) :: radclearcalc(NVIS+1:NVIS+NIR)
96 real(8) :: tbclearcalc(NVIS+1:NVIS+NIR)
97 real(8),allocatable :: radovcalc(:,:)
98 real(8) :: transmsurf(NVIS+1:NVIS+NIR)
99 real(8) :: emiss(NVIS+1:NVIS+NIR)
100 end type avhrr_bgck_iasi
101
102 type(avhrr_bgck_iasi), allocatable :: avhrr_bgck(:) ! avhrr parameters for IASI quality control
103
104contains
105
106
107 !--------------------------------------------------------------------------
108 ! irbg_init
109 !--------------------------------------------------------------------------
110 subroutine irbg_init()
111 !
112 !:Purpose: This subroutine reads the namelist section NAMBGCKIR
113 ! for the module.
114 implicit none
115
116 ! Locals:
117 integer :: nulnam, ierr
118 logical, save :: first = .true.
119 integer, external :: fnom, fclos
120 integer :: instrumentIndex
121 namelist /NAMBGCKIR/ ninst, inst, iwindow, iwindow_alt, ilist1, ilist2, ilist2_pair
122 namelist /NAMBGCKIR/ dtw, dtl, pco2min, pco2max, night_ang, crisCloudFractionThreshold
123
124 if (first) then
125
126 ! set the default values for namelist variables
127 ninst = MPC_missingValue_INT
128 inst(:) = '*EMPTY*'
129 iwindow(:) = 0
130 iwindow_alt(:) = 0
131 ilist1(:,:) = 0
132 ilist2(:,:) = 0
133 ilist2_pair(:,:) = 0
134 dtw = 0.0d0
135 dtl = 0.0d0
136 pco2min = 0.0d0
137 pco2max = 0.0d0
138 night_ang = 0.0d0
139 crisCloudFractionThreshold = -1.d0
140
141 ! read the namelist
142 nulnam = 0
143 ierr = fnom(nulnam, './flnml','FTN+SEQ+R/O', 0)
144 read(nulnam, nml=NAMBGCKIR, iostat=ierr)
145 if (ierr /= 0) call utl_abort('irbg_init: Error reading namelist')
146 if (mmpi_myid == 0) write(*, nml=NAMBGCKIR)
147 if (ninst /= MPC_missingValue_INT) then
148 call utl_abort('irbg_init: check NAMBGCKIR namelist section: ninst should be removed')
149 end if
150 ninst = 0
151 do instrumentIndex = 1, nmaxinst
152 if (inst(instrumentIndex) == '*EMPTY*') exit
153 ninst = ninst + 1
154 end do
155 ierr = fclos(nulnam)
156 first = .false.
157
158 end if
159
160 end subroutine irbg_init
161
162 !--------------------------------------------------------------------------
163 ! irbg_setup
164 !--------------------------------------------------------------------------
165 subroutine irbg_setup()
166 !
167 !:Purpose: Memory allocation for the Hyperspectral Infrared
168 ! background check variables.
169 !
170 implicit none
171
172 ! Locals:
173 integer :: allocStatus(2)
174 integer :: tovsIndex, maxChannelNumber
175 integer :: sensorIndex, channelNumber
176
177 ! Memory allocation for background check related variables
178 allocStatus(:) = 0
179 allocate( tvs_surfaceParameters(tvs_nobtov), stat=allocStatus(1))
180 call utl_checkAllocationStatus(allocStatus(1:1), " irbg_setup tvs_surfaceParameters")
181
182 !___ emissivity by profile
183
184 maxChannelNumber = 1
185 do tovsIndex = 1, tvs_nobtov
186 sensorIndex = tvs_lsensor(tovsIndex)
187 channelNumber = tvs_nchan(sensorIndex)
188 if (channelNumber > maxChannelNumber) maxChannelNumber=channelNumber
189 end do
190
191 allocate( tvs_emissivity (maxChannelNumber,tvs_nobtov), stat=allocStatus(1))
192 call utl_checkAllocationStatus(allocStatus(1:1), " irbg_setup tvs_emissivity")
193
194 do sensorIndex = 1, tvs_nsensors
195 if ( tvs_instruments(sensorIndex) == inst_id_iasi ) then
196 allocate ( avhrr_bgck(tvs_nobtov), stat=allocStatus(1))
197 call utl_checkAllocationStatus(allocStatus(1:1), " irbg_setup avhrr_bgck")
198 exit
199 end if
200 end do
201
202 end subroutine irbg_setup
203
204 !--------------------------------------------------------------------------
205 ! irbg_bgCheckIR
206 !--------------------------------------------------------------------------
207 subroutine irbg_bgCheckIR(columnTrlOnTrlLev, obsSpaceData)
208 !
209 !:Purpose: Do background check on all hyperspectral infrared observations
210 !
211 implicit none
212
213 ! Arguments:
214 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
215 type(struct_obs), intent(inout) :: obsSpaceData
216
217 ! Locals:
218 integer,allocatable :: nobir(:)
219 integer :: headerIndex, idatyp, sensorIndex, instrumentIndex
220 logical :: irDataPresent
221
222 call utl_tmg_start(115,'--BgckInfrared')
223
224 irDataPresent = .false.
225 call obs_set_current_header_list(obsSpaceData,'TO')
226 HEADER0: do
227 headerIndex = obs_getHeaderIndex(obsSpaceData)
228 if (headerIndex < 0) exit HEADER0
229 idatyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
230 if ( tvs_isIdBurpHyperSpectral(idatyp) ) then
231 irDataPresent = .true.
232 end if
233 end do HEADER0
234
235 if ( .not. irDataPresent ) then
236 write(*,*) 'WARNING: WILL NOT RUN irbg_bgCheckIR since no IR'
237 return
238 end if
239
240 write(*,'(A)') " ****************"
241 write(*,'(A)') " BEGIN IR BACKGROUND CHECK"
242 write(*,'(A)') " **************** **************** ****************"
243
244 ! Preliminary initializations
245
246 tvs_nobtov = 0
247 call irbg_init()
248 allocate (nobir(ninst))
249 nobir = 0
250
251
252 ! Loop over all header indices of the 'TO' family
253 ! Set the header list (and start at the beginning of the list)
254 call obs_set_current_header_list(obsSpaceData,'TO')
255 HEADER: do
256 headerIndex = obs_getHeaderIndex(obsSpaceData)
257 if (headerIndex < 0) exit HEADER
258
259 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
260
261 if ( .not. tvs_isIdBurpTovs(idatyp) ) then
262 write(*,*) 'irbg_bgCheckIR: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
263 cycle HEADER ! Proceed to the next header_index
264 end if
265 tvs_nobtov = tvs_nobtov + 1
266 do instrumentIndex=1, ninst
267 if ( tvs_isIdBurpInst(idatyp,inst(instrumentIndex)) ) then
268 nobir(instrumentIndex) = nobir(instrumentIndex) + 1
269 exit
270 end if
271 end do
272 end do HEADER
273
274 do instrumentIndex=1, ninst
275 if (nobir(instrumentIndex) > 0) then
276 do sensorIndex=1,tvs_nsensors
277 if (tvs_instruments(sensorIndex) == tvs_getInstrumentId(inst(instrumentIndex)) ) then
278 call irbg_doQualityControl (columnTrlOnTrlLev, obsSpaceData, inst(instrumentIndex), sensorIndex)
279 end if
280 end do
281 end if
282 end do
283 deallocate (nobir)
284
285 call utl_tmg_stop(115)
286
287 end subroutine irbg_bgCheckIR
288
289 !--------------------------------------------------------------------------
290 ! bgck_get_qcid
291 !--------------------------------------------------------------------------
292 subroutine bgck_get_qcid(instrumentName, qcid )
293 !
294 implicit none
295
296 ! Arguments:
297 character(len=*), intent(in) :: instrumentName
298 integer, intent(out) :: qcid
299
300 ! Locals:
301 integer :: i
302
303 qcid = -1
304
305 do i=1, ninst
306 if (trim(instrumentName) == trim(inst(i))) then
307 qcid = i
308 exit
309 end if
310 end do
311
312 if (qcid == -1) then
313 write(*,*) "Unknown instrument ", instrumentName
314 call utl_abort('bgck_get_qcid')
315 end if
316
317 end subroutine bgck_get_qcid
318
319 !--------------------------------------------------------------------------
320 ! irbg_doQualityControl
321 !--------------------------------------------------------------------------
322 subroutine irbg_doQualityControl (columnTrlOnTrlLev, obsSpaceData, instrumentName, id_opt)
323 !
324 !:Purpose: Quality control of hyperspectral infrared observations.
325 ! assign assimilation flags to observations
326 !
327 implicit none
328
329 ! Arguments:
330 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
331 type(struct_obs), intent(inout) :: obsSpaceData
332 character(len=*), intent(in) :: instrumentName
333 integer, optional, intent(in) :: id_opt
334
335 ! Locals:
336 integer :: jc, nchn, levelIndex, bitIndex, channelNumber, classIndex
337 integer :: columnIndex
338 integer :: iwindo, iwindo_alt
339 integer :: bodyIndex, bodyStart, bodyEnd, headerIndex
340 integer :: idatyp
341 real(8) :: difftop_min
342 integer :: modelTopIndex
343 integer :: count
344 real(8) :: t_effective
345 integer :: allocStatus(27)
346 real(8) :: tg, p0, tskinRetrieved, ptop_T, qs
347 real(8), allocatable :: tt(:), height(:,:)
348 real(8), allocatable :: pressure(:,:)
349 real(8), allocatable :: btObsErr(:), btObs(:), btCalc(:), rcal_clr(:), sfctau(:)
350 real(8), allocatable :: radObs(:), cloudyRadiance(:,:), transm(:,:), emi_sfc(:)
351 real(8), allocatable :: ptop_bt(:), ptop_rd(:)
352 real(8), allocatable :: pmin(:), dtaudp1(:), maxwf(:)
353 real(8), allocatable :: cloudyRadiance_avhrr(:,:)
354 integer, allocatable :: rejflag(:,:)
355 integer, allocatable :: ntop_bt(:), ntop_rd(:)
356 integer, allocatable :: minp(:), fate(:), channelIndexes(:)
357 real(8) :: clfr, sunZenithAngle, satelliteAzimuthAngle, satelliteZenithAngle, sunAzimuthAngle
358 real(8) :: albedo, ice, pcnt_wat, pcnt_reg
359 real(8) :: ptop_eq,ptop_mb
360 real(8) :: ptop_co2(nco2),fcloud_co2(nco2)
361 real(8) :: etop,vtop,ecf,vcf,heff
362 real(8) :: tampon,cfsub
363 real(8) :: tskinRetrieved_avhrr(nClassAVHRR), sfctau_avhrr(NIR), emi_sfc_avhrr(NIR), rcal_clr_avhrr(NIR)
364 real(8) :: ptop_bt_avhrr(NIR,nClassAVHRR), ptop_rd_avhrr(NIR,nClassAVHRR)
365 real(8) :: btObs_avhrr(NIR,nClassAVHRR), radObs_avhrr(NIR,nClassAVHRR), ptop_eq_avhrr(nClassAVHRR)
366 real(8) :: cfrac_avhrr
367 real(8) :: avhrr_surfem1(NIR)
368 real(8) :: albedoThreshold(NIR)
369 integer :: ksurf,ltype
370 integer :: cldflag, lev_start
371 integer :: gncldflag
372 integer :: ichref
373 integer :: ntop_eq, ntop_mb
374 integer :: ngood
375 integer :: ntop_co2(nco2)
376 integer :: cldflag_avhrr(nClassAVHRR),lev_start_avhrr(nClassAVHRR),ichref_avhrr(nClassAVHRR),ntop_rd_avhrr(NIR,nClassAVHRR)
377 integer :: ntop_bt_avhrr(NIR,nClassAVHRR),ntop_eq_avhrr(nClassAVHRR)
378 logical :: assim_all
379 logical :: sunZenithAnglePresent
380 logical :: satelliteAzimuthAnglePresent, satelliteZenithAnglePresent, sunAzimuthAnglePresent
381 integer, parameter :: nn=2
382 integer, parameter :: ilist_avhrr(nn)=(/ 2 ,3 /)
383 integer :: countInvalidChannels
384 logical :: bad
385 real(8), parameter :: sunzenmax=87.12d0
386 real(8) :: minpavhrr(2:3)
387 real(8) :: anisot,zlamb,zcloud,scos,del,deltaphi
388 integer :: ier,ijour,iloc(2:3),co2min(1),co2max(1)
389 integer :: channelIndex,ilist_co2(nco2),ilist_co2_pair(nco2),ilist_he(nch_he)
390 integer :: nlv_T,id,sensorIndex, qcid, nchannels
391 logical :: liasi,lairs,lcris
392 type (rttov_profile), pointer :: profiles(:)
393
394 write (*,*) "Entering irbg_doQualityControl"
395
396 call tvs_getProfile(profiles, 'nl')
397
398 liasi= ( trim(instrumentName) == "IASI" .or. trim(instrumentName) == "iasi")
399 lairs= ( trim(instrumentName) == "AIRS" .or. trim(instrumentName) == "airs")
400 lcris= ( trim(instrumentName) == "CRIS" .or. trim(instrumentName) == "cris" .or. trim(instrumentName) == "CRISFSR" .or. &
401 trim(instrumentName) == "crisfsr" )
402
403 call bgck_get_qcid(instrumentName,qcid)
404
405 if (present(id_opt)) then
406 id = id_opt
407 else
408 ! Find sensor number corresponding to the desired instrument
409 id = -1
410 do sensorIndex = 1, tvs_nsensors
411 if ( trim(tvs_instrumentName(sensorIndex)) == TRIM(instrumentName)) then
412 id = sensorIndex
413 exit
414 end if
415 end do
416 if (id < 0) call utl_abort("irbg_doQualityControl: should not happen !")
417 end if
418
419 ! Find number of profiles
420 count = 0
421
422 ! Loop over all header indices of the 'TO' family
423 call obs_set_current_header_list(obsSpaceData,'TO')
424 HEADER: do
425 headerIndex = obs_getHeaderIndex(obsSpaceData)
426 if (headerIndex < 0) exit HEADER
427 if ( tvs_tovsIndex( headerIndex ) < 0) cycle HEADER
428 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
429 if ( tvs_isIdBurpInst(idatyp,instrumentName) .and. tvs_lsensor(tvs_tovsIndex (headerIndex)) == id ) then
430 count = count + 1
431 end if
432 end do HEADER
433
434 if ( count == 0 ) return
435 ! Find number of channels and RTTOV levels
436
437 nchn = tvs_coefs(id) % coef % fmv_chn
438
439 nlv_T = col_getNumLev(columnTrlOnTrlLev,'TH')
440
441 write(*,*) ' irbg_doQualityControl - nchn ', nchn
442
443
444 ! information to extract (transvidage)
445 !
446 ! tg -- guess skin temperatures (deg K)
447 ! p0 -- surface pressure (hPa)
448 ! tt(nlv_T-1) -- temperature profiles on NWP model levels (deg K)
449 ! height(nlv_T-1,1) -- height profiles on NWP model levels (m)
450 ! qs -- surface specific humidity (kg/kg)
451 ! btObsErr(nchn) -- observation error standard deviation
452 ! btObs(nchn) -- observed brightness temperatures (deg K)
453 ! btCalc(nchn) -- computed brightness temperatures (deg K)
454 ! rcal_clr(nchn) -- computed clear radiances (mw/m2/sr/cm-1)
455 ! sfctau(nchn) -- surface to space transmittances (0-1)
456 ! cloudyRadiance(nchn,nlv_T-1) -- overcast cloudy radiances (mw/m2/sr/cm-1)
457 ! transm(nchn,nlv_T) -- layer to space transmittances (0-1)
458 ! emi_sfc(nchn) -- surface emissivities (0-1)
459 ! ksurf -- surface type in obs file (0, 1)
460 ! clfr -- cloud fraction (%)
461 ! sunZenithAngle -- sun zenith angle (deg)
462 ! satelliteAzimuthAngle -- satellite azimuth angle (deg)
463 ! satelliteZenithAngle -- satellite zenith angle (deg)
464 ! albedo -- surface albedo (0-1)
465 ! ice -- ice fraction (0-1)
466 ! ltype -- surface type (1,...,20)
467 ! pcnt_wat -- water fraction (0-1)
468 ! pcnt_reg -- water fraction in the area (0-1)
469 ! radObs(nchn) -- observed radiances (mW/m2/sr/cm-1)
470
471 allocStatus(:) = 0
472
473 allocate ( btObsErr(nchn), stat= allocStatus(1))
474 allocate ( btObs(nchn), stat= allocStatus(2))
475 allocate ( btCalc(nchn), stat= allocStatus(3))
476 allocate ( rcal_clr(nchn), stat= allocStatus(4))
477 allocate ( sfctau(nchn), stat= allocStatus(5))
478 allocate ( cloudyRadiance(nchn,nlv_T-1), stat= allocStatus(6))
479 allocate ( transm(nchn,nlv_T), stat= allocStatus(7))
480 allocate ( emi_sfc(nchn), stat= allocStatus(8))
481 allocate ( radObs(nchn), stat= allocStatus(11))
482 allocate ( rejflag(nchn,0:bitflag), stat= allocStatus(12))
483 allocate ( ntop_bt(nchn), stat= allocStatus(13))
484 allocate ( ntop_rd(nchn), stat= allocStatus(14))
485 allocate ( ptop_bt(nchn), stat= allocStatus(15))
486 allocate ( ptop_rd(nchn), stat= allocStatus(16))
487 allocate ( minp(nchn), stat= allocStatus(17))
488 allocate ( pmin(nchn), stat= allocStatus(18))
489 allocate ( dtaudp1(nchn), stat= allocStatus(19))
490 allocate ( fate(nchn), stat= allocStatus(20))
491 if (liasi) allocate ( cloudyRadiance_avhrr(NIR,nlv_T-1), stat= allocStatus(21))
492 allocate ( maxwf(nchn), stat= allocStatus(22))
493 allocate ( pressure(nlv_T-1,1), stat= allocStatus(24))
494 allocate ( tt(nlv_T-1), stat= allocStatus(25))
495 allocate ( height(nlv_T-1,1), stat= allocStatus(26))
496 allocate ( channelIndexes(nchn), stat= allocStatus(27))
497 call utl_checkAllocationStatus(allocStatus, " irbg_doQualityControl 1")
498
499 difftop_min = 100000.d0
500 modelTopIndex = 1
501
502 ptop_T = col_getPressure(columnTrlOnTrlLev,1,1,'TH')
503
504 write(*,*) 'TOIT DU MODELE (MB)'
505 write(*,*) 0.01d0 * ptop_T
506 write(*,*) 'NIVEAU DU MODELE DE TRANSFERT RADIATIF LE PLUS PRES DU TOIT DU MODELE'
507 write(*,*) modelTopIndex
508
509 tvs_nobtov = 0
510
511 ! Loop over all header indices of the 'TO' family
512 call obs_set_current_header_list(obsSpaceData, 'TO')
513 HEADER_2: do
514 headerIndex = obs_getHeaderIndex(obsSpaceData)
515 if (headerIndex < 0) exit HEADER_2
516
517 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
518
519 if ( tvs_isIdBurpTovs(idatyp) ) tvs_nobtov = tvs_nobtov + 1
520 if ( tvs_tovsIndex(headerIndex) < 0) cycle HEADER_2
521 if ( tvs_isIdBurpInst(idatyp,instrumentName) .and. tvs_lsensor(tvs_tovsIndex (headerIndex)) == id) then
522 btObs(:) = -1.d0
523 btCalc(:) = -1.d0
524 btObsErr(:) = -1.d0
525 rcal_clr(:) = -1.d0
526 sfctau(:) = -1.d0
527 cloudyRadiance(:,:) = -1.d0
528 transm(:,:) = -1.d0
529 emi_sfc(:) = -1.d0
530 rejflag(:,:) = 0
531 channelIndexes(:) = -1
532
533 if (liasi) then
534 classIndex = 1
535 do columnIndex = OBS_CF1, OBS_CF7
536 avhrr_bgck(headerIndex)%cfrac(classIndex) = obs_headElem_r(obsSpaceData,columnIndex,headerIndex)
537 classIndex = classIndex + 1
538 end do
539 classIndex = 1
540 channelNumber = 1
541 do columnIndex = OBS_M1C1, OBS_M7C6
542 avhrr_bgck(headerIndex)%radmoy(classIndex,channelNumber) = obs_headElem_r(obsSpaceData,columnIndex,headerIndex)
543 channelNumber = channelNumber + 1
544 if (channelNumber > nChanAVHRR) then
545 channelNumber = 1
546 classIndex = classIndex + 1
547 end if
548 end do
549 classIndex = 1
550 channelNumber = 1
551 do columnIndex = OBS_S1C1, OBS_S7C6
552 avhrr_bgck(headerIndex)%radstd(classIndex,channelNumber) = obs_headElem_r(obsSpaceData,columnIndex,headerIndex)
553 channelNumber =channelNumber + 1
554 if (channelNumber> nChanAVHRR) then
555 channelNumber = 1
556 classIndex = classIndex + 1
557 end if
558 end do
559 sunAzimuthAngle = obs_headElem_r(obsSpaceData,OBS_SAZ,headerIndex)
560 sunAzimuthAnglePresent = ( abs(sunAzimuthAngle - MPC_missingValue_R8) > 0.01 )
561 end if
562
563 tg = col_getElem(columnTrlOnTrlLev, 1, headerIndex, 'TG')
564 p0 = col_getElem(columnTrlOnTrlLev, 1, headerIndex, 'P0') * MPC_MBAR_PER_PA_R8
565
566 do levelIndex = 1, nlv_T - 1
567 tt(levelIndex) = col_getElem(columnTrlOnTrlLev, levelIndex+1, headerIndex, 'TT')
568 height(levelIndex,1) = col_getHeight(columnTrlOnTrlLev, levelIndex+1, headerIndex, 'TH')
569 pressure(levelIndex,1)= col_getPressure(columnTrlOnTrlLev, levelIndex+1, headerIndex, 'TH') * MPC_MBAR_PER_PA_R8
570 end do
571 qs = col_getElem(columnTrlOnTrlLev, nlv_T, headerIndex, 'HU')
572
573 bodyStart = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
574 bodyEnd = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) + bodyStart - 1
575 bad = .false.
576 if (lcris) bad=( obs_headElem_i(obsSpaceData, OBS_GQF, headerIndex)/=0 .or. &
577 obs_headElem_i(obsSpaceData, OBS_GQL, headerIndex) /=0)
578 if (liasi) bad=( obs_headElem_i(obsSpaceData, OBS_GQF, headerIndex)/=0 .or. &
579 obs_headElem_i(obsSpaceData, OBS_GQL, headerIndex) >1)
580
581 nchannels = 0 ! number of channels available at that observation point
582 do bodyIndex= bodyStart, bodyEnd
583 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
584 channelNumber = nint(obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex))
585 channelNumber = max( 0,min( channelNumber,tvs_maxChannelNumber + 1 ) )
586 call tvs_getLocalChannelIndexFromChannelNumber(id,channelIndex,channelNumber)
587 if (channelIndex == -1) cycle
588 nchannels = nchannels + 1
589 channelIndexes(nchannels) = channelIndex
590 btObsErr(channelIndex) = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
591 btObs(channelIndex) = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
592
593 ! Flag check on observed BTs ***
594 if (.not.liasi .and. btest(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),2)) rejflag(channelIndex,9) = 1
595 if (bad) rejflag(channelIndex,9) = 1
596
597 ! Gross check on observed BTs ***
598 if (btObs(channelIndex)<150.d0) rejflag(channelIndex,9) = 1
599 if (btObs(channelIndex)>350.d0) rejflag(channelIndex,9) = 1
600 end if
601 end do
602
603 if (nchannels==0) cycle HEADER_2
604
605 do jc = 1, nchannels
606 channelIndex = channelIndexes(jc)
607 if ( channelIndex ==-1) cycle
608 btCalc(channelIndex) = tvs_radiance(tvs_nobtov) % bt(channelIndex)
609 rcal_clr(channelIndex) = tvs_radiance(tvs_nobtov) % clear(channelIndex)
610 sfctau(channelIndex) = tvs_transmission(tvs_nobtov) % tau_total(channelIndex)
611 do levelIndex = 1, nlv_T
612 transm(channelIndex,levelIndex) = tvs_transmission(tvs_nobtov) % tau_levels(levelIndex, channelIndex)
613 end do
614 do levelIndex = 1, nlv_T - 1
615 cloudyRadiance(channelIndex,levelIndex) = tvs_radiance(tvs_nobtov) % overcast(levelIndex, channelIndex)
616 end do
617 emi_sfc(channelIndex) = tvs_emissivity(channelIndex,tvs_nobtov)
618 ! Gross check on computed BTs ***
619 if (btCalc(channelIndex) < 150.d0) rejflag(channelIndex,9) = 1
620 if (btCalc(channelIndex) > 350.d0) rejflag(channelIndex,9) = 1
621 end do
622
623 ksurf = profiles(tvs_nobtov) % skin % surftype
624
625 !Test pour detecter l angle zenithal manquant (-1) ou anormal
626 ! (angle negatif ou superieur a 75 degres )
627 satelliteZenithAngle= obs_headElem_r(obsSpaceData, OBS_SZA, headerIndex)
628 if ( satelliteZenithAngle < 0 .or. satelliteZenithAngle > 75. ) then
629 rejflag(:,9) = 1
630 satelliteZenithAnglePresent = .false.
631 else
632 satelliteZenithAnglePresent = .true.
633 end if
634 clfr = 0.
635 if (lairs .or. lcris) clfr = obs_headElem_r(obsSpaceData, OBS_CLF, headerIndex)
636
637 sunZenithAngle = profiles(tvs_nobtov) % sunzenangle
638 sunZenithAnglePresent = ( abs(sunZenithAngle - MPC_missingValue_R8) > 0.01 )
639
640 if (liasi) then
641 satelliteAzimuthAngle = profiles(tvs_nobtov) % azangle
642 satelliteAzimuthAnglePresent = ( abs(satelliteAzimuthAngle - MPC_missingValue_R8) > 0.01 )
643 end if
644 albedo = tvs_surfaceParameters(tvs_nobtov) % albedo
645 ice = tvs_surfaceParameters(tvs_nobtov) % ice
646 ltype = tvs_surfaceParameters(tvs_nobtov) % ltype
647 if (ltype == 20) ksurf = 2
648 pcnt_wat = tvs_surfaceParameters(tvs_nobtov) % pcnt_wat
649 pcnt_reg = tvs_surfaceParameters(tvs_nobtov) % pcnt_reg
650
651 ! Find TOA radiances converted from observed BT's
652
653 radObs(:) = -1.d0
654
655 channels: do jc = 1, nchannels
656 channelIndex = channelIndexes(jc)
657 if (channelIndex == -1) cycle channels
658 if ( rejflag(channelIndex,9) == 1 ) cycle channels
659 t_effective = tvs_coefs(id) % coef % ff_bco(channelIndex) &
660 + tvs_coefs(id) % coef % ff_bcs(channelIndex) * btObs(channelIndex)
661 radObs(channelIndex) = tvs_coefs(id) % coef % planck1(channelIndex) / &
662 ( exp( tvs_coefs(id) % coef % planck2(channelIndex) / t_effective ) - 1.d0 )
663 end do channels
664
665 ! Set height fields to 'height above ground' fields
666 do levelIndex = 1, nlv_T - 1
667 height(levelIndex,1) = height(levelIndex,1) - height(nlv_T-1,1)
668 end do
669
670 ! ///// ---------------------------------------------------- /////
671 ! ///// Determination of the clear/cloudy profiles (cldflag) /////
672 ! ///// ---------------------------------------------------- /////
673 cldflag = 0
674
675 ! Reference for window channel
676 call tvs_getLocalChannelIndexFromChannelNumber(id, iwindo, iwindow(qcid) )
677 call tvs_getLocalChannelIndexFromChannelNumber(id, iwindo_alt, iwindow_alt(qcid) )
678 ichref = -1
679 if (iwindo /= -1) then
680 if ( rejflag(iwindo,9) == 0) then
681 ichref = iwindo
682 end if
683 end if
684 if (ichref == -1 .and. iwindo_alt /= -1) then
685 if ( rejflag(iwindo_alt,9) == 0) then
686 ichref = iwindo_alt
687 end if
688 end if
689
690 if (ichref == -1) then
691 cldflag = -1
692 rejflag(:,9) = 1
693 write(*,*) 'WARNING'
694 write(*,*) 'WINDOW AND ALTERNATE WINDOW CHANNEL OBSERVATIONS'
695 write(*,*) 'HAVE BEEN REJECTED. '
696 write(*,*) 'ALL '//instrumentName//' OBSERVATIONS FROM THIS PROFILE REJECTED'
697 end if
698
699 ! -- Cloud top based on matching observed brightness temperature
700 ! -- at a reference surface channel with background temperature profile (ptop_eq)
701 ! -- on guess vertical levels.
702
703 lev_start = 0
704
705 !iopt2=1 : calcul de la hauteur en hPa ptop_mb et du ntop_mb correspondant
706 call cloud_height ( ptop_mb, ntop_mb, btObs, cldflag, tt, &
707 height(:,1), p0, pressure(:,1), ichref, lev_start, iopt2 )
708
709 !iopt1=2 : calcul de la hauteur em metres ptop_eq et du ntop_eq correspondant
710 call cloud_height ( ptop_eq, ntop_eq, btObs, cldflag, tt, &
711 height(:,1), p0, pressure(:,1), ichref, lev_start, iopt1 )
712
713 if (liasi) then
714 ! appel de RTTOV pour calculer les radiances des 3 canaux IR (3b, 4 et 5) de AVHRR 3
715
716 call get_avhrr_emiss(emi_sfc(channelIndexes(1:nchannels)),tvs_coefs(id) % coef % ff_cwn(channelIndexes(1:nchannels)), &
717 nchannels,avhrr_surfem1)
718
719 call tovs_rttov_avhrr_for_IASI(headerIndex,avhrr_surfem1,tvs_satellites(id))
720
721 !The value computed will be used only if sunZenithAnglePresent is true
722 call convert_avhrr(sunZenithAngle, avhrr_bgck(headerIndex) )
723 call stat_avhrr(avhrr_bgck(headerIndex))
724
725 lev_start_avhrr(:) = 0
726 cldflag_avhrr(:) = 0
727 do classIndex = 1, nClassAVHRR
728 btObs_avhrr(:,classIndex) = avhrr_bgck(headerIndex) % tbmoy(classIndex,:)
729 radObs_avhrr(1:NIR,classIndex) = avhrr_bgck(headerIndex) % radmoy(classIndex,NVIS + 1:NIR + NVIS)
730 rcal_clr_avhrr(:) = avhrr_bgck(headerIndex) % radclearcalc(:)
731 emi_sfc_avhrr(:) = avhrr_bgck(headerIndex) % emiss(:)
732 sfctau_avhrr(:) = avhrr_bgck(headerIndex) % transmsurf(:)
733
734 do levelIndex = 1, nlv_T - 1
735 cloudyRadiance_avhrr(:,levelIndex) = avhrr_bgck(headerIndex) % radovcalc(levelIndex,:)
736 end do
737
738 if (btObs_avhrr(2,classIndex) > 100.d0 ) then
739 ichref_avhrr(classIndex) = 2
740 else if (btObs_avhrr(3,classIndex) > 100.d0 ) then
741 ichref_avhrr(classIndex) = 3
742 else
743 ichref_avhrr(classIndex) = -1
744 cldflag_avhrr(classIndex) = -1
745 end if
746
747 call cloud_height (ptop_eq_avhrr(classIndex),ntop_eq_avhrr(classIndex), btObs_avhrr(:,classIndex),cldflag_avhrr(classIndex),tt, &
748 height(:,1),p0,pressure,ichref_avhrr(classIndex),lev_start_avhrr(classIndex),iopt1)
749 end do
750
751 end if
752
753 ! -- Clear/cloudy profile detection using the garand & nadon algorithm
754
755 call garand1998nadon (cldflag, btObs,tg,tt, &
756 height(:,1),ptop_eq,ntop_eq,ichref)
757
758 if (liasi) then
759 do classIndex=1,nClassAVHRR
760 call garand1998nadon (cldflag_avhrr(classIndex), btObs_avhrr(:,classIndex),tg,tt, &
761 height(:,1),ptop_eq_avhrr(classIndex),ntop_eq_avhrr(classIndex),ichref_avhrr(classIndex))
762 end do
763 end if
764
765 ! Further tests to remove potential cloudy profiles
766 ! Test # A
767 ! In daytime, set cloudy if cloud fraction over 5%
768 cfsub = -1.d0
769 if (lairs) then
770 if ( cldflag == 0 .and. clfr > 5.d0 .and. sunZenithAngle < 90.d0 .and. sunZenithAnglePresent) then
771 cldflag = 1
772 cfsub = 0.01d0 * clfr !conversion % -> 0-1
773 end if
774 end if
775 if (lcris) then
776 if ( cldflag == 0 .and. clfr >= 0.d0 .and. clfr > crisCloudFractionThreshold &
777 .and. crisCloudFractionThreshold >= 0.d0 ) then
778 cldflag = 1
779 cfsub = 0.01d0 * clfr !conversion % -> 0-1
780 end if
781 end if
782 ! Test # B
783 ! Set cloudy if temperature difference between guess (tg)
784 ! and estimated true (tskinRetrieved) skin temperatures is over threshold
785
786 call estim_ts(tskinRetrieved, tg, emi_sfc, rcal_clr, radObs, &
787 sfctau, cldflag, ichref, tvs_coefs(id) )
788
789 if ( cldflag == 0 .and. ksurf == 1 &
790 .and. abs(tskinRetrieved-tg) > dtw ) cldflag = 1
791
792 if ( cldflag == 0 .and. ksurf /= 1 &
793 .and. abs(tskinRetrieved-tg) > dtl ) cldflag = 1
794
795 if (liasi) then
796
797 do classIndex = 1, nClassAVHRR
798 call estim_ts(tskinRetrieved_avhrr(classIndex), tg,emi_sfc_avhrr,rcal_clr_avhrr,radObs_avhrr(:,classIndex), &
799 sfctau_avhrr,cldflag_avhrr(classIndex),ichref_avhrr(classIndex), coefs_avhrr)
800 end do
801
802 do classIndex = 1, nClassAVHRR
803 if ( cldflag_avhrr(classIndex) == 0 .and. ksurf == 1 &
804 .and. abs(tskinRetrieved_avhrr(classIndex)-tg) > dtw ) cldflag_avhrr(classIndex) = 1
805
806 if ( cldflag_avhrr(classIndex) == 0 .and. ksurf /= 1 &
807 .and. abs(tskinRetrieved_avhrr(classIndex)-tg) > dtl ) cldflag_avhrr(classIndex) = 1
808
809 end do
810
811 !criteres AVHRR utilisant les canaux visibles (de jour seulement)
812 if ( sunZenithAngle < sunzenmax .and. sunZenithAnglePresent .and. satelliteAzimuthAnglePresent .and. &
813 sunAzimuthAnglePresent .and. satelliteZenithAnglePresent) then
814 anisot = 1.d0
815
816 if (albedo < albedoThresholdQC) then
817 deltaphi = abs(satelliteAzimuthAngle - sunAzimuthAngle )
818 if (deltaphi > 180.d0) deltaphi = 360.d0 - deltaphi
819 call visocn(sunZenithAngle,satelliteZenithAngle,deltaphi,anisot,zlamb,zcloud,IER)
820 albedoThreshold = 10.d0 * max(1.d0,anisot)
821 else
822 albedoThreshold = 100.d0 * albedo + 10.d0
823 end if
824
825 if (anisot < 1.5d0) then !to avoid sun glint
826 scos = cos ( sunZenithAngle * MPC_DEGREES_PER_RADIAN_R8 )
827 call cor_albedo (del, scos )
828 albedoThreshold = albedoThreshold * del
829 do classIndex = 1, nClassAVHRR
830 if (avhrr_bgck(headerIndex)%albedmoy(classIndex,1) > albedoThreshold(1) ) then
831 cldflag_avhrr(classIndex) = 1
832 end if
833 !static AVHRR thresholds v3
834 do channelIndex = 1, NVIS
835 if (avhrr_bgck(headerIndex)%albedmoy(classIndex,channelIndex) > seuilalb_static(channelIndex,ksurf) ) then
836 cldflag_avhrr(classIndex) = 1
837 end if
838 end do
839 end do
840
841 end if
842 end if
843
844 ! Calcul de la pseudo fraction nuageuse AVHRR
845
846 cfrac_avhrr = 0.d0
847 do classIndex = 1, nClassAVHRR
848 if (cldflag_avhrr(classIndex) == 1) cfrac_avhrr = cfrac_avhrr + avhrr_bgck(headerIndex) % cfrac(classIndex)
849 end do
850
851 cfsub = -1.0d0
852 if ( cldflag == 0 .and. cfrac_avhrr > 5.d0 ) then
853 cldflag = 1
854 cfsub = 0.01d0 * min(cfrac_avhrr, 100.d0) !conversion % -> 0-1 avec seuil car parfois cfrac_avhrr=101
855 end if
856
857 !AVHRR Homogeneity criteria
858 if (cldflag == 0 .and. sunZenithAnglePresent) then
859 ijour = 1
860 if (sunZenithAngle < 90.d0) ijour=2
861 ! 1 NUIT
862 ! 2 JOUR
863 if (ijour == 2) then
864 do channelIndex=1,NVIS
865 if (avhrr_bgck(headerIndex)%albstd_pixeliasi(channelIndex) > seuilalb_homog(channelIndex,ksurf) ) cldflag = 1
866 end do
867 end if
868 do channelIndex=NVIS+1,NVIS+NIR
869 if (avhrr_bgck(headerIndex)%tbstd_pixelIASI(channelIndex) > seuilbt_homog(channelIndex,ksurf,ijour)) cldflag = 1
870 end do
871 end if
872 end if
873
874 gncldflag = cldflag
875
876 ! ///// ------------------------------------------------------- /////
877 ! ///// DETERMINATION OF THE ASSIMILABLE OBSERVATIONS (rejflag) /////
878 ! ///// ------------------------------------------------------- /////
879
880
881 ! -- FIRST TestS TO REJECT OBSERVATIONS
882
883
884 ! *** Test # 1 ***
885 ! *** Do not assimilate where cloudy ***
886
887 if ( cldflag == 1 ) then
888 rejflag(:,11) = 1
889 rejflag(:,23) = 1
890 end if
891
892 ! *** Test # 2 ***
893 ! *** Gross check on valid BTs ***
894
895 ! already done
896
897
898 ! -- Cloud top based on matching
899 ! -- observed brightness temperature with background temperature profiles (ptop_bt)
900 ! -- or computed observed radiances with background radiance profiles (ptop_rd)
901 ! -- on rttov vertical levels
902
903 lev_start = 0
904
905 do channelIndex = 1, nch_he
906 call tvs_getLocalChannelIndexFromChannelNumber(id,ilist_HE(channelIndex),ilist1(qcid,channelIndex))
907 end do
908 !here the case for which one of the channelindex is -1 is treated in cloud_top
909 call cloud_top ( ptop_bt,ptop_rd,ntop_bt,ntop_rd, &
910 btObs,tt,height(:,1),rcal_clr,p0,radObs,cloudyRadiance,pressure(:,1), &
911 cldflag, lev_start, iopt2, ihgt, ilist_he,rejflag_opt=rejflag,ichref_opt=ichref)
912
913 if (liasi) then
914 lev_start_avhrr(:) = 0
915 do classIndex = 1, nClassAVHRR
916 call cloud_top( ptop_bt_avhrr(:,classIndex),ptop_rd_avhrr(:,classIndex),ntop_bt_avhrr(:,classIndex),ntop_rd_avhrr(:,classIndex), &
917 btObs_avhrr(:,classIndex),tt,height(:,1),rcal_clr_avhrr,p0,radObs_avhrr(:,classIndex),cloudyRadiance_avhrr,pressure(:,1), &
918 cldflag_avhrr(classIndex),lev_start_avhrr(classIndex),iopt2,ihgt,ilist_avhrr)
919 end do
920 end if
921
922 ! -- reference channel for co2-slicing
923
924 do channelIndex = 1, nco2
925 call tvs_getLocalChannelIndexFromChannelNumber(id, ilist_co2(channelIndex), ilist2(qcid,channelIndex) )
926 call tvs_getLocalChannelIndexFromChannelNumber(id, ilist_co2_pair(channelIndex), ilist2_pair(qcid,channelIndex) )
927 end do
928
929 countInvalidChannels = 0
930 do channelIndex = 1, nco2
931 if ( ilist_co2(channelIndex) == -1 .or. ilist_co2_pair(channelIndex) == -1 ) then
932 countInvalidChannels = countInvalidChannels + 1
933 else if ( rejflag(ilist_co2(channelIndex),9) == 1 .or. &
934 rejflag(ilist_co2_pair(channelIndex),9) == 1 ) then
935 countInvalidChannels = countInvalidChannels + 1
936 end if
937 end do
938
939 if (countInvalidChannels == nco2) then
940 cldflag = -1
941 rejflag(:,9) = 1
942 write(*,*) 'WARNING'
943 write(*,*) 'CO2 REFERENCE AND ALTERNATE CHANNEL OBSERVATIONS'
944 write(*,*) 'HAVE BEEN REJECTED. '
945 write(*,*) 'ALL '//instrumentName//' OBSERVATIONS FROM THIS PROFILE REJECTED'
946 end if
947
948 ! Equivalent height of selected window channel
949 heff = MPC_missingValue_R8
950 call tvs_getLocalChannelIndexFromChannelNumber(id,channelIndex,ilist1(qcid,2))
951 if (channelIndex /= -1) heff = ptop_rd( channelIndex )
952
953 if (ichref == iwindo_alt) then
954 call tvs_getLocalChannelIndexFromChannelNumber(id,channelIndex,ilist1(qcid,3))
955 if (channelIndex /= -1) heff = ptop_rd( channelIndex )
956 end if
957 ! Cloud top based on co2 slicing
958
959 co2min = minloc( abs( pressure(:,1) - pco2min ) )
960 co2max = minloc( abs( pressure(:,1) - pco2max ) )
961
962 lev_start = max( min(lev_start,co2max(1)), co2min(1) )
963
964 call co2_slicing ( ptop_co2, ntop_co2, fcloud_co2, &
965 rcal_clr, cloudyRadiance, radObs, p0, pressure(:,1), cldflag, rejflag, &
966 lev_start, ichref, ilist_co2, ilist_co2_pair)
967
968 ! -- Find consensus cloud top and fraction
969
970 call seltop ( etop,vtop,ecf,vcf,ngood, heff,ptop_co2,fcloud_co2, &
971 cfsub,ptop_mb,p0,cldflag,gncldflag )
972
973 if (liasi) then
974 ! Correction pour les nuages trop bas:
975 ! en principe Pco2 < Heff.
976 ! on cherche les cas pathologiques avec Pco2>Min(Heff(AVHRR))
977 minpavhrr(2:3) = 12200
978 iloc(2:3) = -1 ! pour eviter les catastrophes...
979 do classIndex=1,nClassAVHRR
980 if (avhrr_bgck(headerIndex)%cfrac(classIndex) > 0.d0) then
981 if (ptop_rd_avhrr(2,classIndex) < minpavhrr(2)) then
982 iloc(2) = classIndex
983 minpavhrr(2) = ptop_rd_avhrr(2,classIndex)
984 end if
985 if (ptop_rd_avhrr(3,classIndex) < minpavhrr(3)) then
986 iloc(3) = classIndex
987 minpavhrr(3) = ptop_rd_avhrr(3,classIndex)
988 end if
989 end if
990 end do
991 if ( iloc(2) /= -1 .and. iloc(3) /= -1) then ! pour eviter les catastrophes...
992 ! on se limite aux cas "surs" ou les deux hauteurs effectives sont > a Pco2
993 ! et ou un accord raisonnable existe entre les deux hauteurs effectives
994 if ( iloc(2) == iloc(3) .and. &
995 minpavhrr(2) < etop .and. &
996 minpavhrr(3) < etop .and. &
997 abs(minpavhrr(2)- minpavhrr(3)) < 25.d0 .and. &
998 cldflag_avhrr(iloc(2)) /= -1 .and. cldflag_avhrr(iloc(3)) /= -1) then
999
1000 if (ecf == 0.d0 .and. cldflag == 1) then
1001 ! cas predetermine nuageux mais ramene a clair
1002 ecf = 0.01d0 * min(100.d0,cfrac_avhrr)
1003 ! cette ligne peut generer des fractions nuageuses inferieures a 20 %.
1004 etop = 0.5d0 * (minpavhrr(2) + minpavhrr(3))
1005 end if
1006
1007 if (ecf > 0.d0 .and. cldflag == 1) then
1008 !cas predetermine nuageux pas ramene clair (==normal)
1009 etop = 0.5d0 * ( minpavhrr(2) + minpavhrr(3))
1010 end if
1011
1012 if (cldflag == 0) then
1013 !cas predetermine clair ... que faire
1014 cldflag = 1
1015 etop = 0.5d0 * (minpavhrr(2) + minpavhrr(3))
1016 ecf = 0.01d0 * min(100.d0,cfrac_avhrr)
1017 end if
1018 end if
1019 end if
1020 end if
1021
1022 ! -- Find minimum level of sensitivity for channel assimilation not sensible to clouds
1023 call min_pres_new (maxwf, minp, pmin, dtaudp1, p0, transm(:,2:nlv_T), pressure(:,1), cldflag, modelTopIndex)
1024 ! -- ASSIMILATION OF OBSERVATIONS WHEN CLOUDY PROFILES
1025
1026 ! *** Test # 3 ***
1027 ! *** Assimilation above clouds (refinement of test 1) ***
1028 ! *** Set security margin to 2x the std on height from CO2-slicing ***
1029
1030 tampon = max(50.d0, 2.d0*vtop)
1031
1032 do channelIndex = 1, nchn
1033 if ( rejflag(channelIndex,11) == 1 .and. rejflag(channelIndex,23) == 1 .and. etop - tampon > pmin(channelIndex) ) then
1034 rejflag(channelIndex,11) = 0
1035 rejflag(channelIndex,23) = 0
1036 end if
1037 end do
1038
1039 ! Look at the fate of the observations
1040 fate(:) = sum(rejflag(:,:), DIM=2)
1041
1042 ! Further reasons to reject observations
1043 do channelIndex = 1, nchn
1044
1045 if ( fate(channelIndex) == 0 ) then
1046
1047 ! *** Test # 4 ***
1048 ! *** Background check, do not assimilate if O-P > 3sigma ***
1049
1050 if ( abs(btObs(channelIndex) - btCalc(channelIndex)) > 3.d0 * btObsErr(channelIndex) ) then
1051 rejflag(channelIndex,9) = 1
1052 rejflag(channelIndex,16) = 1
1053 end if
1054
1055 ! *** Test # 5 ***
1056 ! *** Do not assimilate shortwave channels during the day ***
1057
1058 if ( tvs_coefs(id) % coef % ff_cwn(channelIndex) >= 2000.d0 .and. sunZenithAngle < night_ang .and. sunZenithAnglePresent) then
1059 rejflag(channelIndex,11) = 1
1060 rejflag(channelIndex,7) = 1
1061 end if
1062
1063 ! *** Test # 6 ***
1064 ! *** Do not assimilate surface channels over land ***
1065
1066 if ( minp(channelIndex) == nlv_T .or. p0-pmin(channelIndex) < 100.d0 ) then
1067 if ( ksurf == 0 ) then
1068 rejflag(channelIndex,11) = 1 !!! comment this line if assimilation under conditions
1069 rejflag(channelIndex,19) = 1 !!! comment this line if assimilation under conditions
1070 if ( pcnt_wat > 0.01d0 .or. pcnt_reg > 0.1d0 .or. emi_sfc(channelIndex) < 0.97d0 ) then
1071 rejflag(channelIndex,11) = 1
1072 rejflag(channelIndex,19) = 1
1073 end if
1074
1075 ! *** Test # 7 ***
1076 ! *** Do not assimilate surface channels over water under conditions ***
1077
1078 else if ( ksurf == 1 ) then
1079 if ( pcnt_wat < 0.99d0 .or. pcnt_reg < 0.97d0 .or. &
1080 ice > 0.001d0 .or. albedo >= albedoThresholdQC .or. emi_sfc(channelIndex) < 0.9d0 ) then
1081 rejflag(channelIndex,11) = 1
1082 rejflag(channelIndex,19) = 1
1083 end if
1084
1085 ! *** Test # 8 ***
1086 ! *** Do not assimilate surface channels over sea ice ***
1087
1088 else if ( ksurf == 2 ) then
1089 rejflag(channelIndex,11) = 1
1090 rejflag(channelIndex,19) = 1
1091 end if
1092 end if
1093
1094 end if
1095
1096 ! *** Test # 9 ***
1097 ! *** Do not assimilate if jacobian has a significant contribution over model top ***
1098
1099 ! Condition valid if model top at 10mb or lower only
1100 if ( nint(ptop_T) >= 1000 ) then
1101 if ( rejflag(channelIndex,9) /= 1 .and. dtaudp1(channelIndex) > 0.50d0 ) then
1102 rejflag(channelIndex,11) = 1
1103 rejflag(channelIndex,21) = 1
1104 end if
1105 end if
1106
1107 ! Condition valid if model top at 10mb or lower only
1108 if ( nint(ptop_T) >= 1000 ) then
1109 if ( rejflag(channelIndex,9) /= 1 .and. transm(channelIndex,1) < 0.99d0 ) then
1110 rejflag(channelIndex,11) = 1
1111 rejflag(channelIndex,21) = 1
1112 end if
1113 end if
1114
1115 ! Condition valid if model top is higher than 10 mb
1116 ! with computation made on model levels instead of RTTOV levels
1117 ! this test is a litle bit more strict than before
1118 ! decrease of the order of 1-2% in the number of high peaking radiance IR assimilated
1119 ! not dramatic so could be kept as is
1120 if ( nint(ptop_T) < 1000 ) then
1121 if ( rejflag(channelIndex,9) /= 1 .and. transm(channelIndex,1) < 0.95d0 ) then
1122 rejflag(channelIndex,11) = 1
1123 rejflag(channelIndex,21) = 1
1124 end if
1125 end if
1126
1127 end do
1128
1129 nchannels =0
1130 do bodyIndex= bodyStart, bodyEnd
1131 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
1132 nchannels = nchannels + 1
1133 if (btest(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),8)) rejflag(channelIndexes(nchannels),8) = 1
1134 end if
1135 end do
1136
1137 ! For each profile, are all non-blacklisted channels assimilated
1138
1139 assim_all = .true.
1140 fate(:) = sum(rejflag(:,:),DIM=2)
1141
1142 chn: do channelIndex = 1, nchn
1143 if ( rejflag(channelIndex,8) == 0 ) then
1144 if ( fate(channelIndex) /= 0 ) then
1145 assim_all = .false.
1146 exit chn
1147 end if
1148 end if
1149 end do chn
1150
1151 if (.not.assim_all) then
1152 call obs_headSet_i(obsSpaceData, OBS_ST1, headerIndex,ibset(obs_headElem_i(obsSpaceData,OBS_ST1,headerIndex),6) )
1153 end if
1154 ! -- Addition of background check parameters to burp file
1155
1156 call obs_headSet_r(obsSpaceData, OBS_ETOP, headerIndex, etop )
1157 call obs_headSet_r(obsSpaceData, OBS_VTOP, headerIndex, vtop )
1158 call obs_headSet_r(obsSpaceData, OBS_ECF, headerIndex, 100.d0 * ecf )
1159 call obs_headSet_r(obsSpaceData, OBS_VCF, headerIndex, 100.d0 * vcf )
1160 call obs_headSet_r(obsSpaceData, OBS_HE, headerIndex, heff )
1161 call obs_headSet_r(obsSpaceData, OBS_ZTSR, headerIndex, tskinRetrieved )
1162 call obs_headSet_i(obsSpaceData, OBS_NCO2, headerIndex, ngood)
1163 call obs_headSet_r(obsSpaceData, OBS_ZTM, headerIndex, tt(nlv_T-1) )
1164 call obs_headSet_r(obsSpaceData, OBS_ZTGM, headerIndex, tg )
1165 call obs_headSet_r(obsSpaceData, OBS_ZLQM, headerIndex, qs )
1166 call obs_headSet_r(obsSpaceData, OBS_ZPS, headerIndex, 100.d0 * p0 )
1167 call obs_headSet_i(obsSpaceData, OBS_STYP, headerIndex, ksurf )
1168
1169 do bodyIndex = bodyStart, bodyEnd
1170 channelNumber = nint(obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex))
1171 channelNumber = max(0, min(channelNumber, tvs_maxChannelNumber + 1))
1172 call tvs_getLocalChannelIndexFromChannelNumber(id,channelIndex,channelNumber)
1173 if (channelIndex == -1) cycle
1174 call obs_bodySet_r(obsSpaceData,OBS_SEM,bodyIndex,emi_sfc(channelIndex))
1175 do bitIndex = 0, bitflag
1176 if ( rejflag(channelIndex,bitIndex) == 1 ) &
1177 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex,IBSET(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),bitIndex))
1178 end do
1179 end do
1180
1181 end if
1182
1183 end do HEADER_2
1184
1185 deallocate ( channelIndexes, stat= allocStatus(1))
1186 deallocate ( height, stat= allocStatus(2))
1187 deallocate ( tt, stat= allocStatus(3))
1188 deallocate ( pressure, stat= allocStatus(4))
1189 deallocate ( maxwf, stat= allocStatus(5))
1190 if (liasi) deallocate ( cloudyRadiance_avhrr , stat= allocStatus(6))
1191 deallocate ( fate, stat= allocStatus(7))
1192 deallocate ( dtaudp1, stat= allocStatus(8))
1193 deallocate ( pmin, stat= allocStatus(9))
1194 deallocate ( minp, stat= allocStatus(10))
1195 deallocate ( ptop_rd, stat= allocStatus(11))
1196 deallocate ( ptop_bt, stat= allocStatus(12))
1197 deallocate ( ntop_rd, stat= allocStatus(13))
1198 deallocate ( ntop_bt, stat= allocStatus(14))
1199 deallocate ( rejflag, stat= allocStatus(15))
1200 deallocate ( radObs, stat= allocStatus(16))
1201 deallocate ( emi_sfc, stat= allocStatus(17))
1202 deallocate ( transm, stat= allocStatus(18))
1203 deallocate ( cloudyRadiance, stat= allocStatus(19))
1204 deallocate ( sfctau, stat= allocStatus(20))
1205 deallocate ( rcal_clr, stat= allocStatus(21))
1206 deallocate ( btCalc, stat= allocStatus(22))
1207 deallocate ( btObs, stat= allocStatus(23))
1208 deallocate ( btObsErr, stat= allocStatus(24))
1209 call utl_checkAllocationStatus(allocStatus, " irbg_doQualityControl", .false.)
1210 nullify(profiles)
1211
1212 end subroutine irbg_doQualityControl
1213
1214
1215 !--------------------------------------------------------------------------
1216 ! convert_avhrr
1217 !--------------------------------------------------------------------------
1218 subroutine convert_avhrr(sunzen, avhrr)
1219 !
1220 !:Purpose: conversion des radiance IR en temperatures de brillance
1221 ! et des radiances visibles en "albedo"
1222
1223 implicit none
1224
1225 ! Arguments:
1226 real(8), intent(in) :: sunzen ! Solar zenith angle
1227 type (avhrr_bgck_iasi), intent(inout) :: avhrr ! Structure containing AVHRR observations
1228
1229 ! Locals:
1230 integer :: classIndex
1231 real(8) :: bt(NIR), dbtsdrad(NIR)
1232 real(8) :: freq(NIR), offset(NIR), slope(NIR)
1233
1234 freq = coefs_avhrr%coef%ff_cwn (:)
1235 offset = coefs_avhrr%coef%ff_bco(:)
1236 slope = coefs_avhrr%coef%ff_bcs(:)
1237
1238 do classIndex=1, nClassAVHRR
1239 call calcbt(avhrr % radmoy(classIndex,4:6), bt, dbtsdrad, freq, offset, slope)
1240 avhrr % tbmoy(classIndex,4:6) = bt(1:3)
1241 avhrr % tbstd(classIndex,4:6) = avhrr % radstd(classIndex,4:6) * dbtsdrad(1:3)
1242 call calcreflect(avhrr % radmoy(classIndex,1:3), sunzen, avhrr % albedmoy(classIndex,1:3) )
1243 call calcreflect(avhrr % radstd(classIndex,1:3), sunzen, avhrr % albedstd(classIndex,1:3) )
1244 end do
1245
1246 end subroutine convert_avhrr
1247
1248 !--------------------------------------------------------------------------
1249 ! calcreflect
1250 !--------------------------------------------------------------------------
1251 subroutine calcreflect(rad, sunzen, reflect)
1252 !
1253 !:Purpose: Computes Top of Atmosphere Albedo as defined by equation (4)
1254 ! of Rao et al. Int. J. of Remote Sensing, 2003, vol 24, no 9, 1913-1924
1255 !
1256 implicit none
1257
1258 ! Arguments:
1259 real(8), intent(in) :: rad(nvis) ! radiances array
1260 real(8), intent(in) :: sunzen ! Sun zenith angle
1261 real(8), intent(out):: reflect(nvis) ! TOA albedo en %
1262
1263 ! Locals:
1264 real(8),parameter :: solar_filtered_irradiance(nvis) = (/139.873215d0,232.919556d0,14.016470d0/)
1265 ! equivalent widths, integrated solar irradiance, effective central wavelength
1266 !0.084877,139.873215,0.632815
1267 !0.229421,232.919556,0.841679
1268 !0.056998,14.016470,1.606119
1269 real(8) :: radb ! radiance en W/m2/str
1270 integer :: i
1271
1272 do i = 1, nvis
1273 if (rad(i) >= 0.0d0 ) then
1274 radb = rad(i) / 1000.0d0
1275 reflect(i) = (MPC_PI_R8 * radb) / solar_filtered_irradiance(i)
1276 if (sunzen < 90.0d0 ) reflect(i) = reflect(i) / cos(sunzen * MPC_RADIANS_PER_DEGREE_R8)
1277 else
1278 reflect(i) = -1
1279 end if
1280 end do
1281
1282 end subroutine calcreflect
1283
1284 !--------------------------------------------------------------------------
1285 ! calcbt
1286 !--------------------------------------------------------------------------
1287 subroutine calcbt(rad, tb, dtbsdrad, freq, offset, slope)
1288 !
1289 !:Purpose: Computes brightness temperature (bt) and the first derivative of
1290 ! bt with respect to radiance, from radiance, frequencies
1291 !
1292 !
1293 implicit none
1294
1295 ! Arguments:
1296 real(8), intent(in) :: rad(:) ! Radiance
1297 real(8), intent(in) :: freq(size(rad)) ! Channel wavenumber (cm-1)
1298 real(8), intent(in) :: offset(size(rad)) !
1299 real(8), intent(in) :: slope(size(rad)) !
1300 real(8), intent(out) :: tb(size(rad)) ! Brightness Temperature
1301 real(8), intent(out) :: dtbsdrad(size(rad)) ! Derivative of tb wrt radiance
1302
1303 ! Locals:
1304 integer :: channelIndex, nchan
1305 real(8) :: radtotal, tstore, planck1, planck2
1306 real(8), parameter :: c1= 1.19106590d-05 ! First Planck constant
1307 real(8), parameter :: c2= 1.438833d0 ! Second Planck constant
1308
1309 nchan = size(rad)
1310
1311 do channelIndex = 1, nchan
1312 if (rad(channelIndex) > 1.d-20) then
1313 planck2 = c2 * freq(channelIndex)
1314 planck1 = c1 * ( freq(channelIndex) ** 3 )
1315 tstore = planck2 / log( 1.0d0 + planck1 / rad(channelIndex) )
1316 tb(channelIndex) = ( tstore - offset(channelIndex) ) / slope(channelIndex)
1317
1318 radtotal = rad(channelIndex)
1319
1320 dtbsdrad(channelIndex) = planck1 * tstore ** 2 / ( planck2 * radtotal * ( radtotal + planck1 ) )
1321
1322 dtbsdrad(channelIndex) = dtbsdrad(channelIndex) / slope(channelIndex)
1323
1324 else
1325 tb(channelIndex) = 0.d0
1326 dtbsdrad(channelIndex) = 0.d0
1327 end if
1328
1329 end do
1330
1331 end subroutine calcbt
1332
1333 !--------------------------------------------------------------------------
1334 ! stat_avhrr
1335 !--------------------------------------------------------------------------
1336 subroutine stat_avhrr( avhrr )
1337 !
1338 !:Purpose: calcul de statistiques
1339 ! sur l'information sous-pixel AVHRR
1340 !
1341 implicit none
1342
1343 ! Arguments:
1344 type (avhrr_bgck_iasi), intent(inout) :: avhrr
1345
1346 ! Locals:
1347 integer :: classIndex, channelIndex
1348 real(8) :: sumFrac(nvis+nir), sumBt(nvis+1:nvis+nir),sumBt2(nvis+1:nvis+nir)
1349 real(8) :: sumAlb(1:nvis),sumAlb2(1:nvis)
1350
1351 sumFrac(:) = 0.d0
1352 sumBt(:) = 0.d0
1353 sumBt2(:) = 0.d0
1354 sumAlb(:) = 0.d0
1355 sumAlb2(:) = 0.d0
1356
1357 do classIndex = 1, nClassAVHRR
1358 if (avhrr%cfrac(classIndex) > 0.d0 ) then
1359 do channelIndex = 1, NVIS
1360 if (avhrr%albedmoy(classIndex,channelIndex) >= 0.d0 ) then
1361 sumFrac(channelIndex) = sumFrac(channelIndex) + avhrr%cfrac(classIndex)
1362 sumAlb(channelIndex) = sumAlb(channelIndex) + avhrr%cfrac(classIndex) * avhrr%albedmoy(classIndex,channelIndex)
1363 sumAlb2(channelIndex) = sumAlb2(channelIndex) + avhrr%cfrac(classIndex) * ( avhrr%albedmoy(classIndex,channelIndex)**2 + avhrr%albedstd(classIndex,channelIndex)**2)
1364 end if
1365 end do
1366 do channelIndex = 1+NVIS, NVIS+NIR
1367 if (avhrr%tbmoy(classIndex,channelIndex) > 0.d0 ) then
1368 sumFrac(channelIndex) = sumFrac(channelIndex) + avhrr%cfrac(classIndex)
1369 sumBt(channelIndex) = sumBt(channelIndex) + avhrr%cfrac(classIndex) * avhrr%tbmoy(classIndex,channelIndex)
1370 sumBt2(channelIndex) = sumBt2(channelIndex) + avhrr%cfrac(classIndex) * (avhrr%tbmoy(classIndex,channelIndex)**2 + avhrr%tbstd(classIndex,channelIndex)**2 )
1371 end if
1372 end do
1373 end if
1374 end do
1375
1376 do channelIndex = 1, NVIS
1377 if (sumFrac(channelIndex) > 0.d0 ) then
1378 sumAlb(channelIndex) = sumAlb(channelIndex) / sumFrac(channelIndex)
1379 sumAlb2(channelIndex) = sumAlb2(channelIndex)/sumFrac(channelIndex) - sumAlb(channelIndex)**2
1380 if (sumAlb2(channelIndex) > 0.d0) then
1381 sumAlb2(channelIndex) = sqrt( sumAlb2(channelIndex) )
1382 else
1383 sumAlb2(channelIndex) = 0.d0
1384 end if
1385 end if
1386 end do
1387
1388 do channelIndex = NVIS+1, NVIS+NIR
1389 if (sumFrac(channelIndex) > 0.d0 ) then
1390 sumBt(channelIndex) = sumBt(channelIndex) / sumFrac(channelIndex)
1391 sumBt2(channelIndex) = sumBt2(channelIndex)/sumFrac(channelIndex) - sumBt(channelIndex)**2
1392 if (sumBt2(channelIndex) > 0.d0) then
1393 sumBt2(channelIndex)= sqrt ( sumBt2(channelIndex) )
1394 else
1395 sumBt2(channelIndex) = 0.d0
1396 end if
1397 end if
1398 end do
1399
1400 avhrr % tbstd_pixelIASI = sumBt2
1401 avhrr % albstd_pixeliasi = sumAlb2
1402
1403 end subroutine stat_avhrr
1404
1405 !--------------------------------------------------------------------------
1406 ! co2_slicing
1407 !--------------------------------------------------------------------------
1408 subroutine co2_slicing (ptop, ntop, fcloud, &
1409 rcal, cloudyRadiance, radObs, p0, plev, cldflag, rejflag, &
1410 lev_start, ichref, ilist, ilist_pair)
1411 !
1412 !:Purpose: cloud top height computation.
1413 ! cloud top from co2 slicing and cloud fraction estimate
1414 !
1415 implicit none
1416
1417 ! Arguments:
1418 real(8), intent(in) :: rcal(:) ! computed clear radiances (mW/m2/sr/cm-1)
1419 real(8), intent(in) :: plev(:) ! pressure levels (hPa)
1420 real(8), intent(in) :: cloudyRadiance(size(rcal),size(plev)) ! computed cloud radiances from each level (mW/m2/sr/cm-1)
1421 real(8), intent(in) :: radObs(size(rcal)) ! Observed radiances (mW/m2/sr/cm-1)
1422 real(8), intent(in) :: p0 ! surface pressure (hPa)
1423 integer, intent(in) :: ichref ! window channel to predetermine clear
1424 integer, intent(in) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
1425 integer, intent(in) :: rejflag(1:,0:) ! flags for rejected observations
1426 integer, intent(in) :: ilist(nco2) ! first list of channels
1427 integer, intent(in) :: ilist_pair(nco2) ! second list channe
1428 integer, intent(inout) :: lev_start ! Level to start iteration (ideally tropopause)
1429 real(8), intent(out) :: ptop(nco2) ! Cloud top (hPa)
1430 real(8), intent(out) :: fcloud(nco2) ! Cloud fraction
1431 integer, intent(out) :: ntop(nco2) ! Nearest pressure level corresponding to ptop (ptop <= p0)
1432
1433 ! Locals
1434 integer :: j,jch,jc,jpmax,jmax, nlev,nchn
1435 integer :: sumrej
1436 real(8) :: rapg
1437 real(8),allocatable :: drap(:,:), a_drap(:), fc(:,:)
1438 real(8) :: val,val1,val2,val3,fcint
1439 real(8) :: emi_ratio
1440 integer :: jc_pair
1441 integer :: iter,niter
1442 real(8), parameter :: eps = 1.D-12
1443
1444 ptop(:) = -1.d0
1445 ntop(:) = -1
1446 fcloud(:) = -1.d0
1447
1448 ! Profile not assimilated if data from 2 windows channels bad
1449 ! and/or if data from 2 reference co2 channels bad
1450
1451 if ( cldflag == -1 ) return
1452
1453 nlev = size(plev)
1454 nchn = size(rcal)
1455
1456 ! Define closest level jpmax to surface pressure p0
1457
1458 jpmax = nlev
1459
1460 do J = lev_start, nlev
1461 if ( plev(J) > p0 ) then
1462 jpmax = J
1463 exit
1464 end if
1465 end do
1466
1467 ! define jmax as last level for co2-slicing calculations
1468
1469 jmax = jpmax - 1
1470
1471 ! predetermined clear window channel, all nco2 estimates clear
1472
1473 sumrej = sum(rejflag(ichref,:))
1474
1475 if ( sumrej == 0 ) then
1476 ptop(:) = p0
1477 ntop(:) = jpmax
1478 fcloud(:) = 0.d0
1479 return
1480 end if
1481
1482 allocate(fc(nchn,nlev), drap(nco2,nlev), a_drap(nlev) )
1483
1484 channels: do jch = 1, nco2
1485
1486 jc = ilist(jch)
1487 jc_pair = ilist_pair(jch)
1488 fc(jc_pair,:) = rcal(jc_pair) - cloudyRadiance(jc_pair,:)
1489 niter = 1
1490 if ( jch > 13) niter = 2
1491
1492 iteration: do iter = 1, niter
1493 drap(jch,:) = 9999.d0
1494 ntop(jch) = -1
1495 ! calcul emi_ratio
1496 if (jch > 13) then
1497 if ( iter == 1 ) then
1498 emi_ratio = 1.0376d0
1499 else
1500 emi_ratio = 1.09961d0 - 0.09082d0 * fcloud(jch)
1501 end if
1502 else
1503 emi_ratio = 1.0d0
1504 end if
1505
1506 fc(jc,:) = rcal(jc) - cloudyRadiance(jc,:)
1507
1508 ! Gross check failure
1509
1510 if ( rejflag(jc,9) == 1 ) cycle channels
1511 if ( rejflag(jc_pair,9) == 1 ) cycle channels
1512
1513 if ( abs( rcal(jc_pair) - radObs(jc_pair) ) > eps ) then
1514 rapg = (rcal(jc) - radObs(jc)) / (rcal(jc_pair) - radObs(jc_pair))
1515 else
1516 rapg = 0.0d0
1517 end if
1518
1519 do J = lev_start, jpmax
1520 if ( fc(jc,J) > 0.d0 .and. fc(jc_pair,J) > 0.d0 ) &
1521 drap(jch,J) = rapg - (fc(jc,J) / fc(jc_pair,J)) * emi_ratio
1522 end do
1523
1524 a_drap(:) = abs( drap(jch,:) )
1525
1526 levels: do J = lev_start + 1, jmax
1527
1528 ! * do not allow fc negative (i.e. drap(jch,j) = 9999.)
1529
1530 if ( drap(jch,J) > 9000.d0 .and. &
1531 a_drap(J-1) < eps .and. &
1532 a_drap(J+1) < eps ) cycle channels
1533
1534 val = drap(jch,J) / drap(jch,J - 1)
1535
1536 ! * find first, hopefully unique, zero crossing
1537
1538 if ( val < 0.d0 ) then
1539
1540 ! * conditions near zero crossing of isolated minimum need monotonically
1541 ! * decreasing drap from j-3 to j-1 as well increasing from j to j+1
1542
1543 val1 = drap(jch,J - 2) / drap(jch,J - 1)
1544 val2 = drap(jch,J - 3) / drap(jch,J - 1)
1545 val3 = drap(jch,J) / drap(jch,J + 1)
1546
1547 if ( val1 > 0.d0 .and. &
1548 val2 > 0.d0 .and. &
1549 val3 > 0.d0 .and. &
1550 a_drap(J-2) > a_drap(J-1) .and. &
1551 a_drap(J-3) > a_drap(J-2) .and. &
1552 a_drap(J) < 9000.d0 .and. &
1553 a_drap(J+1) > a_drap(J) ) &
1554 then
1555 ptop(jch) = plev(J)
1556 ntop(jch) = J
1557 end if
1558
1559 exit levels
1560
1561 end if
1562
1563 end do levels
1564
1565 j = ntop(jch)
1566
1567 ! * special cases of no determination
1568
1569 if ( j < 1) then
1570 ptop(jch) = -1.d0
1571 ntop(jch) = -1
1572 fcloud(jch) = -1.d0
1573 cycle channels
1574 end if
1575
1576 if ( J <= lev_start .or. drap(jch,J) > 9000.d0 ) then
1577 !if ( iter == 1) then
1578 ptop(jch) = -1.d0
1579 ntop(jch) = -1
1580 fcloud(jch) = -1.d0
1581 !end if
1582 cycle channels
1583 end if
1584
1585 if ( abs( cloudyRadiance(jc,J) - rcal(jc) ) > 0.d0 ) &
1586 fcloud(jch) = (radObs(jc) - rcal(jc)) / &
1587 (cloudyRadiance(jc,J) - rcal(jc))
1588
1589 ! * find passage to zero if it exists and interpolate to exact pressure
1590
1591 ptop(jch) = plev(J - 1) - drap(jch,J - 1) / &
1592 ( drap(jch,J) - drap(jch,J-1) ) * ( plev(J) - plev(J-1) )
1593 ! * find cloud radiance at zero crossing to use to get cloud fraction
1594
1595 fcint = fc(jc,J - 1) + ( fc(jc,J) - fc(jc,J - 1) ) / &
1596 ( plev(J) - plev(J - 1) ) * ( ptop(jch) - plev(J - 1) )
1597
1598 ! * find cloud fraction based on exact cloud top
1599
1600 if ( abs(fcint) > 0.d0 ) &
1601 fcloud(jch) = ( rcal(jc) - radObs(jc) ) / fcint
1602
1603 fcloud(jch) = min ( fcloud(jch), 1.5d0 )
1604 fcloud(jch) = max ( fcloud(jch), -0.5d0 )
1605
1606 if (fcloud(jch) < 0.0d0 .or. fcloud(jch) > 1.0d0 ) cycle channels
1607
1608 end do iteration
1609
1610 end do channels
1611
1612 deallocate(fc, drap, a_drap )
1613
1614 end subroutine co2_slicing
1615
1616 !--------------------------------------------------------------------------
1617 ! seltop
1618 !--------------------------------------------------------------------------
1619 subroutine seltop (etop, vtop, ecf, vcf, ngood, he, ht, cf, cfsub, ptop_mb, p0, cldflag, gncldflag)
1620 !
1621 !:Purpose: Select cloud top by averaging co2-slicing results
1622 ! judged correct. all missing values are -1.
1623 !
1624 implicit none
1625
1626 ! Arguments:
1627 real(8), intent(in) :: he ! equivalent cloud top heights from a window channel (hPa)
1628 real(8), intent(in) :: ht(nco2) ! cloud tops from co2-slicing (hPa)
1629 real(8), intent(in) :: cf(nco2) ! effective cloud fraction for co2-slicing
1630 real(8), intent(in) :: p0 ! surface pressure in (hPa)
1631 real(8), intent(in) :: cfsub ! visible ("subpixel") cloud fraction
1632 integer, intent(in) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
1633 integer, intent(in) :: gncldflag! Garand Nadon cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
1634 real(8), intent(out) :: etop ! consensus cloud top (hPa)
1635 real(8), intent(out) :: vtop ! corresponding variance on etop (hPa)
1636 real(8), intent(out) :: ecf ! consensus effective cloud fraction
1637 real(8), intent(out) :: vcf ! corresponding variance on ecf
1638 integer, intent(out) :: ngood ! number of good estimates
1639 real(8), intent(in) :: ptop_mb ! height (mb) from cloud_height subroutine
1640
1641 ! Locals:
1642 integer :: n,jch
1643 real(8) :: H(nco2),F(nco2)
1644
1645 etop = -1.d0
1646 vtop = -1.d0
1647 ecf = -1.d0
1648 vcf = -1.d0
1649 ngood = 0
1650
1651 ! Profile not assimilated if data from 2 windows channels bad
1652 ! and/or if data from 2 reference co2 channels bad
1653 if ( cldflag == -1 ) return
1654
1655 n = 0
1656 H(:) = 0.d0
1657 F(:) = 0.d0
1658
1659 do jch = 1, nco2
1660
1661 ! Check for zero cloud fraction
1662
1663 if ( CF(jch) > -0.9d0 .and. CF(jch) < 1.D-6 ) then
1664 n = n + 1
1665 H(n) = p0
1666 F(n) = 0.d0
1667 else
1668
1669
1670 ! Consider only valid values of cloud fraction above some threshold
1671
1672 ! Important logic: for values above 1.0 of co2-slicing cloud fraction,
1673 ! set it to 1.0 and force the top equal to the effective height he.
1674 ! co2-slicing not allowed to give estimates below he, which happens
1675 ! for cloud fraction cf > 1.0.
1676
1677 if ( ht(jch) > 0.0d0 ) then
1678 n = n + 1
1679 H(n) = ht(jch)
1680 F(n) = min(CF(jch), 1.0d0)
1681 F(n) = max(F(n), 0.d0)
1682 if ( CF(jch) > 1.0d0 ) H(n) = he
1683 end if
1684 end if
1685
1686 end do
1687
1688
1689 ngood = n
1690
1691 ! Compute mean and variance
1692
1693 if ( n >= 1 ) then
1694
1695 call calcul_median_fast(n,H,F,etop,ecf)
1696
1697 vtop = sqrt ( sum((H(1:n) - etop)**2) / n )
1698 vcf = sqrt ( sum((F(1:n) - ecf)**2) / n )
1699
1700 if ( n == 1 ) then
1701 vtop = 50.d0
1702 vcf = 0.20d0
1703 end if
1704
1705 else
1706
1707 ! If no solution from co2-slicing, and not predetermined clear,
1708 ! assume cloudy with top equal to effective height he;
1709 ! however if he is very close to surface pressure p0, assume clear.
1710
1711 etop = he
1712 ecf = 1.0d0
1713 if (cfsub >= 0.05d0) then
1714 ecf = cfsub
1715 etop = min( min(he,ptop_mb) , p0 - 50.0d0)
1716 end if
1717 vtop = 50.d0
1718 vcf = 0.30d0
1719 if ( he > (p0 - 10.d0) ) ecf = 0.d0
1720 if ( gncldflag == 0 ) then
1721 ecf = 0.0d0
1722 etop = p0
1723 end if
1724 end if
1725
1726 if ( ecf < 0.05d0 ) then
1727 ecf = 0.0d0
1728 etop = p0
1729 end if
1730
1731 end subroutine seltop
1732
1733 !--------------------------------------------------------------------------
1734 ! calcul_median_fast
1735 !--------------------------------------------------------------------------
1736 subroutine calcul_median_fast(nEstimates, Hin, Fin, ctp, cfr)
1737 !
1738 !:Purpose: Compute cloud fraction and height median.
1739 !
1740 !
1741 implicit none
1742
1743 ! Arguments:
1744 integer, intent(in) :: nEstimates ! Number of Co2 slicing estimates
1745 real(8), intent(in) :: Hin(:) ! Array of Height estimates
1746 real(8), intent(in) :: Fin(:) ! Array of Cloud fraction estimates
1747 real(8), intent(out) :: ctp ! Median of cloud top pressure estimates
1748 real(8), intent(out) :: cfr ! Corresponding cloud fraction
1749
1750 ! Locals:
1751 integer :: index(nEstimates)
1752 real(4) :: H(nEstimates)
1753 integer :: i
1754
1755 if (nEstimates == 1) then
1756 ctp = Hin(nEstimates)
1757 cfr = Fin(nEstimates)
1758 else
1759 H(1:nEstimates) = Hin(1:nEstimates)
1760 call IPSORT(index,H,nEstimates)
1761 if (mod(nEstimates,2) == 0) then ! N - pair
1762 i = index(nEstimates / 2)
1763 else ! N - impair
1764 i = index(1 + nEstimates / 2)
1765 end if
1766 ctp = Hin(i)
1767 cfr = Fin(i)
1768 end if
1769
1770 end subroutine calcul_median_fast
1771
1772 !--------------------------------------------------------------------------
1773 ! min_pres_new
1774 !--------------------------------------------------------------------------
1775 subroutine min_pres_new( maxheight, minp, pmin, dt1, p0, tau, plev, cldflag, &
1776 modelTopIndex )
1777 !
1778 !:Purpose: from total transmittance array, find minimum height
1779 ! level of sensitivity for a number of profiles and channels.
1780 ! this may be used to select for assimilation only the
1781 ! observations without sensitivity to clouds, that is the
1782 ! response function significant only above cloud level.
1783 ! the criterion is that dtau/dplev > 0.01 for a 100 mb layer.
1784 !
1785 !
1786 implicit none
1787
1788 ! Arguments:
1789 real(8), intent(out) :: maxHeight(:) ! Height (hPa) of the maximum of the weighting function
1790 integer, intent(out) :: minp(size(maxHeight)) ! vertical level corresponding to pmin
1791 real(8), intent(out) :: pmin(size(maxHeight)) ! minimum height of sensitivity (hPa)
1792 real(8), intent(out) :: dt1(size(maxHeight)) ! value of 'dtau/dlogp' at model top
1793 real(8), intent(in) :: p0 ! surface pressure (hPa)
1794 real(8), intent(in) :: plev(:) ! pressure levels (hPa)
1795 real(8), intent(in) :: tau(size(maxHeight),size(plev)) ! layer to space transmittances (0.-1.)
1796 integer, intent(in) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
1797 integer, intent(in) :: modelTopIndex ! rt model level nearest to model top
1798
1799 ! Locals:
1800 real(8) :: maxwf
1801 integer :: levelIndex, channelIndex, ipos(1), nlev, nchn
1802 real(8),allocatable :: wfunc(:), rap(:)
1803
1804 minp(:) = -1
1805 pmin(:) = -1.d0
1806 dt1(:) = -1.d0
1807
1808 if ( cldflag == -1 ) return
1809
1810 nlev = size(plev)
1811 nchn = size(maxHeight)
1812
1813 allocate( wfunc(nlev-1), rap(nlev-1) )
1814
1815 do levelIndex = 1, nlev - 1
1816 rap(levelIndex) = log( plev(levelIndex + 1) / plev(levelIndex) )
1817 end do
1818
1819 channels: do channelIndex = 1, nchn
1820
1821 ! Profile not assimilated if data from 2 windows channels bad
1822 ! and/or if data from 2 reference co2 channels bad
1823
1824 do levelIndex = 1, nlev
1825 if ( tau(channelIndex,levelIndex) < 0.d0) cycle channels
1826 end do
1827
1828 minp(channelIndex) = nlev
1829 pmin(channelIndex) = min(plev(nlev),p0)
1830
1831 ! Compute entire array of dtau/dlog(P)
1832
1833 do levelIndex = 1, nlev - 1
1834 wfunc(levelIndex) = (tau(channelIndex,levelIndex) - tau(channelIndex,levelIndex + 1)) / rap(levelIndex)
1835 end do
1836
1837 dt1(channelIndex) = wfunc(modelTopIndex)
1838
1839 ! If channel sees the surface, don't recalculate minp and pmin
1840
1841 if ( tau(channelIndex,nlev) > 0.01d0 ) cycle channels
1842
1843 ! Recherche du maximum
1844 ipos = maxloc( wfunc(:) )
1845 ! Calcul de la valeur du maximum
1846 maxwf = wfunc(ipos(1))
1847 ! maximum entre les 2 niveaux puisque WF calculee pour une couche finie ( discutable ?)
1848 maxheight(channelIndex)= 0.5d0 * ( plev(ipos(1)) + plev(ipos(1) + 1) )
1849
1850 ! If channel doesn't see the surface, see where dtau/dlog(plev) becomes important
1851 ! for recomputation of minp and pmin.
1852
1853 do levelIndex = nlev - 1, ipos(1), -1
1854 if ( ( wfunc(levelIndex)/ maxwf ) > 0.01d0) then
1855 minp(channelIndex) = levelIndex + 1
1856 pmin(channelIndex) = min(plev(levelIndex + 1),p0)
1857 exit
1858 end if
1859 end do
1860
1861 end do channels
1862
1863 deallocate(wfunc, rap )
1864
1865 end subroutine min_pres_new
1866
1867 !--------------------------------------------------------------------------
1868 ! cloud_height
1869 !--------------------------------------------------------------------------
1870 subroutine cloud_height ( ptop, ntop, btObs, cldflag, tt, height, p0, plev, &
1871 ichref, lev_start, iopt )
1872 !
1873 !:Purpose:
1874 ! Computation of cloud top height (above the ground)
1875 ! based on matching observed brightness temperature at a
1876 ! reference surface channel with background temperature profile.
1877 ! to use with one reference channel. used here on model levels.
1878 !
1879 implicit none
1880
1881 ! Arguments:
1882 real(8), intent(out) :: ptop ! Chosen equivalent cloud tops (in hpa|m with iopt = 1|2)
1883 integer, intent(out) :: ntop ! Number of possible ptop solutions
1884 real(8), intent(in) :: btObs(:) ! Observed brightness temperature (deg k)
1885 integer, intent(in) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
1886 real(8), intent(in) :: tt(:) ! Temperature profiles (deg K)
1887 real(8), intent(in) :: height(size(tt))! Height profiles above ground (m)
1888 real(8), intent(in) :: p0 ! Surface pressure (hPa)
1889 real(8), intent(in) :: plev(size(tt)) ! Pressure levels (hPa)
1890 integer, intent(in) :: ichref ! Chosen reference surface channel
1891 integer, intent(inout) :: lev_start ! Level to start iteration (ideally tropopause)
1892 integer, intent(in) :: iopt ! Levels using plev (1) or height (2)
1893
1894 ! Locals:
1895 integer :: itop
1896 integer :: nht, nlev
1897 real(8), allocatable :: ht(:)
1898
1899 nlev = size(tt)
1900 allocate(ht(nlev))
1901
1902 if ( iopt == 1 ) then
1903
1904 ptop = p0
1905 ntop = 1
1906
1907 if ( cldflag == -1 ) return
1908
1909 call get_top ( ht, nht, btObs(ichref), tt, plev, lev_start, iopt )
1910
1911 itop = 1
1912 if ( nht >= 2 ) itop = 2
1913 ptop = min ( ht(itop), p0 )
1914 ntop = nht
1915
1916 else if ( iopt == 2 ) then
1917
1918 ptop = 0.d0
1919 ntop = 1
1920
1921 if ( cldflag == -1 ) return
1922
1923 call get_top ( ht,nht, btObs(ichref),tt,height,lev_start,iopt )
1924
1925 itop = 1
1926 if ( nht >= 2 ) itop = 2
1927 ptop = max ( ht(itop), 0.d0 )
1928 ntop = nht
1929
1930 end if
1931
1932 deallocate(ht)
1933
1934 end subroutine cloud_height
1935
1936 !--------------------------------------------------------------------------
1937 ! garand1998nadon
1938 !--------------------------------------------------------------------------
1939 subroutine garand1998nadon ( cldflag, btObs, tg, tt, height, &
1940 ptop_eq, ntop_eq, ichref )
1941 !
1942 !:Purpose: Determine if the profiles are clear or cloudy based on
1943 ! the algorithm of garand & nadon 98 j.clim v11 pp.1976-1996
1944 ! with channel iref.
1945 implicit none
1946
1947 ! Arguments:
1948 integer, intent(inout) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
1949 real(8), intent(in) :: btObs(:) ! Observed brightness temperatures (K)
1950 real(8), intent(in) :: tg ! Guess skin temperatures (K)
1951 real(8), intent(in) :: tt(:) ! Guess temperature profiles (K)
1952 real(8), intent(in) :: height(size(tt))! Guess height profile above ground (m)
1953 real(8), intent(in) :: ptop_eq ! Chosen equivalent cloud tops (m)
1954 integer, intent(in) :: ntop_eq ! Number of possible ptop_eq solutions
1955 integer, intent(in) :: ichref ! Chosen reference surface channel
1956
1957 ! Locals:
1958 integer :: ninv
1959 real(8) :: lev(2)
1960
1961 lev(1) = 222.d0
1962 lev(2) = 428.d0
1963
1964 if ( cldflag == -1 ) return
1965
1966 if ( btObs(ichref) >= tg - 3.d0 .and. btObs(ichref) <= tg + 3.d0 ) then
1967 cldflag = 0
1968 return
1969 end if
1970
1971 if ( btObs(ichref) >= tg - 4.d0 .and. btObs(ichref) <= tg - 3.d0 ) then
1972 if ( ptop_eq > 1100.d0 ) then
1973 cldflag = 1
1974 return
1975 else
1976 cldflag = 0
1977 return
1978 end if
1979 end if
1980
1981 if ( ptop_eq > 728.d0 ) then
1982 cldflag = 1
1983 return
1984 end if
1985
1986 if ( tg - btObs(ichref) > 8.d0 ) then
1987 if ( ntop_eq >= 3 ) then
1988 if ( ptop_eq > 73.d0 ) then
1989 cldflag=1
1990 return
1991 else
1992 cldflag=0
1993 return
1994 end if
1995 else
1996 call monotonic_inversion (ninv, tg,tt,height,lev(1))
1997 if ( ninv == 1 ) then
1998 if ( ptop_eq > 222.d0 ) then
1999 cldflag = 1
2000 return
2001 else
2002 cldflag = 0
2003 return
2004 end if
2005 else
2006 cldflag = 0
2007 return
2008 end if
2009 end if
2010 end if
2011
2012 if ( tg - btObs(ichref) > 5.d0 ) then
2013 if ( ntop_eq >= 3 ) then
2014 if ( ptop_eq > 222.d0 ) then
2015 cldflag = 1
2016 return
2017 else
2018 cldflag = 0
2019 return
2020 end if
2021 else
2022 call monotonic_inversion (ninv, tg,tt,height,lev(2))
2023 if ( ninv == 1) then
2024 if( ptop_eq > 428.d0 ) then
2025 cldflag = 1
2026 return
2027 else
2028 cldflag = 0
2029 return
2030 end if
2031 else
2032 cldflag = 0
2033 end if
2034 end if
2035 else
2036 cldflag = 0
2037 end if
2038
2039 end subroutine garand1998nadon
2040
2041 !--------------------------------------------------------------------------
2042 ! monotonic_inversion
2043 !--------------------------------------------------------------------------
2044 subroutine monotonic_inversion ( ninvr, tg, tt, height, lvl )
2045 !
2046 !:Purpose: Determine if there is a presence (ninvr=1) or not (ninvr=0)
2047 ! of a temperature inversion going from the surface up to the
2048 ! height lvl
2049 !
2050 !:Arguments:
2051 ! :ninvr: PRESENCE (1) OR NOT (0) OF A TEMPERATURE INVERSION
2052 ! FROM THE SURFACE TO HEIGHT LVL
2053 !
2054 implicit none
2055
2056 ! Arguments:
2057 real(8), intent(in) :: tt(:) ! Temperature profile (K)
2058 real(8), intent(in) :: height(size(tt))! Height profile above ground (m)
2059 real(8), intent(in) :: tg ! Skin temperature (K)
2060 real(8), intent(in) :: lvl ! Height to search for temperature inversion (m)
2061 integer, intent(out) :: ninvr ! Number of inversions
2062
2063 ! Locals:
2064 integer :: levelIndex, nlevels
2065
2066 ninvr = 0
2067 nlevels = size ( tt )
2068 if ( tg - tt(nlevels) < 0.d0 ) then
2069 ninvr = 1
2070 do levelIndex = nlevels - 1, 1, -1
2071 if ( height(levelIndex) > lvl ) exit
2072 if ( tt(levelIndex+1) - tt(levelIndex) > 0.d0 ) then
2073 ninvr = 0
2074 exit
2075 end if
2076 end do
2077 end if
2078
2079 end subroutine monotonic_inversion
2080
2081
2082 !--------------------------------------------------------------------------
2083 ! estim_ts
2084 !--------------------------------------------------------------------------
2085 subroutine estim_ts( ts, tg, emi, rcal, radobs, sfctau, cldflag, &
2086 ichref, myCoefs )
2087 !
2088 !:Purpose: Get an estimated skin temperature by inversion of
2089 ! radiative transfer equation assuming guess t and q profiles
2090 ! are perfect. designed for a single channel ichref and nprf
2091 ! profiles. assumes a real tg (guess) over oceans and a tg
2092 ! with hypothesis of unity emissivity over land.
2093 !
2094 !:Note: Uses rcal = B(TG)*EMI*SFCtau + ATMOS_PART
2095 ! ts = B(ts)*EMI*SFCtau + ATMOS_PART
2096 ! SOLVES FOR ts
2097 !
2098 implicit none
2099
2100 ! Arguments:
2101 integer, intent(in) :: ichref ! Reference surface channel (subset values)
2102 integer, intent(in) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
2103 real(8), intent(in) :: tg ! Guess skin temperature (K)
2104 real(8), intent(in) :: emi(:) ! Surface emissivities from window channel (0.-1.)
2105 real(8), intent(in) :: rcal(size(emi)) ! Computed clear radiances (mW/m2/sr/cm-1)
2106 real(8), intent(in) :: radobs(size(emi)) ! Observed radiances (mW/m2/sr/cm-1)
2107 real(8), intent(in) :: sfctau(size(emi)) ! Surface to space transmittances (0.-1.)
2108 real(8), intent(out) :: ts ! Retrieved skin temperature (-1. for missing)
2109 type(rttov_coefs), intent(in) :: myCoefs ! RTTOV coefficients structure
2110
2111 ! Locals:
2112 real(8) :: rtg,radtg
2113 real(8) :: radts,tstore,t_effective
2114
2115 ts = -1.d0
2116
2117 if ( cldflag /= 0 ) return
2118 if ( ichref == -1 ) return
2119
2120
2121 ! Transform guess skin temperature to plank radiances
2122
2123 t_effective = myCoefs % coef % ff_bco(ichref) + myCoefs % coef % ff_bcs(ichref) * TG
2124
2125 radtg = myCoefs % coef % planck1(ichref) / &
2126 ( exp( myCoefs % coef % planck2(ichref) / t_effective ) - 1.0d0 )
2127
2128
2129 ! Compute TOA planck radiances due to guess skin planck radiances
2130
2131 rtg = radtg * EMI(ichref) * SFCtau(ichref)
2132
2133
2134 ! Compute true skin planck radiances due to TOA true planck radiances
2135
2136 radts = ( RADOBS(ichref) + rtg - rcal(ichref) ) / &
2137 ( EMI(ichref) * SFCtau(ichref) )
2138
2139 if (radts <= 0.d0) then
2140 write(*,'(A25,1x,8e14.6)') "Warning, negative radts", RADOBS(ichref), rtg, rcal(ichref), EMI(ichref), SFCtau(ichref), &
2141 ( RADOBS(ichref) + rtg - rcal(ichref) ), ( EMI(ichref) * SFCtau(ichref) ), radts
2142 write(*,*) "Skipping tskin retrieval."
2143 return
2144 end if
2145
2146
2147 ! Transform true skin planck radiances to true skin temperatures
2148
2149 tstore = myCoefs % coef % planck2(ichref) / log( 1.0d0 + myCoefs % coef % planck1(ichref) / radts )
2150
2151 ts = ( tstore - myCoefs % coef % ff_bco(ichref) ) / myCoefs % coef % ff_bcs(ichref)
2152
2153
2154 end subroutine estim_ts
2155
2156 !--------------------------------------------------------------------------
2157 ! cloud_top
2158 !--------------------------------------------------------------------------
2159 subroutine cloud_top ( ptop_bt, ptop_rd, ntop_bt, ntop_rd, btObs, &
2160 tt, height, rcal, p0, radObs, cloudyRadiance, plev, &
2161 cldflag, lev_start, iopt, ihgt, &
2162 ilist,rejflag_opt, ichref_opt )
2163 !
2164 !:Purpose: Computation of cloud top height (above the ground)
2165 ! based on matching observed brightness temperature with
2166 ! background temperature profiles and/or computed observed
2167 ! radiances with background radiance profiles.
2168 ! to use with more than one channel. used here on rttov levels.
2169 !
2170 !
2171 implicit none
2172
2173 ! Arguments:
2174 real(8), intent(out) :: ptop_bt(:) ! Chosen equivalent cloud tops based on brightness temperatures (in hpa|m with iopt = 1|2)
2175 real(8), intent(out) :: ptop_rd(:) ! Chosen equivalent cloud tops based on radiances (in hpa|m with iopt = 1|2)
2176 integer, intent(out) :: ntop_bt(:) ! Number of possible ptop_bt solutions
2177 integer, intent(out) :: ntop_rd(:) ! Number of possible ptop_rd solutions
2178 real(8), intent(in) :: btObs(:) ! Observed brightness temperautres (K)
2179 real(8), intent(in) :: tt(:) ! Temperature profiles (K)
2180 real(8), intent(in) :: height(:) ! Height profiles above ground (m)
2181 real(8), intent(in) :: rcal(:) ! Computed clear radiances (mW/m2/sr/cm-1)
2182 real(8), intent(in) :: p0 ! Surface pressure (hPa)
2183 real(8), intent(in) :: radObs(:) ! Computed observed radiances (mW/m2/sr/cm-1)
2184 real(8), intent(in) :: cloudyRadiance(:,:) ! Computed cloud radiances from each level (hPa)
2185 real(8), intent(in) :: plev(:) ! Pressure levels (hPa)
2186 integer, intent(in) :: cldflag ! Cloudy flag (0 Clear, 1 Cloudy, -1 undefined)
2187 integer, intent(inout) :: lev_start ! Level to start iteration (ideally tropopause)
2188 integer, intent(in) :: iopt ! Levels using plev (1) or height (2)
2189 integer, intent(in) :: ihgt ! Get *_bt* only (0), *_rd* only (1), both (2)
2190 integer, intent(in) :: ilist(:) ! List of the channel numbers (subset values)
2191 integer, optional, intent(in) :: rejflag_opt(1:,0:)! Flags for rejected observations
2192 integer, optional, intent(in) :: ichref_opt ! Reference surface channel (subset value)
2193
2194 ! Locals:
2195 integer :: jch,jc,itop,nht,i10,i,nlev,nch
2196 real(8),allocatable :: ht(:)
2197 logical :: clear, cloudy
2198
2199 ! Profile not assimilated if data from 2 windows channels bad
2200
2201 ptop_bt(:) = -10.d0
2202 ptop_rd(:) = -10.d0
2203 ntop_bt(:) = 0.d0
2204 ntop_rd(:) = 0.d0
2205
2206 if ( cldflag == -1 ) return
2207
2208 nlev = size( tt )
2209 i10=1
2210 do i=2, nlev
2211 if (plev(i - 1) <= 100.d0 .and. plev(i) > 100.d0) then
2212 i10 = i
2213 exit
2214 end if
2215 end do
2216
2217 ! predetermined clear
2218 if ( present(rejflag_opt)) then
2219 clear = ( sum( rejflag_opt(ichref_opt,:) ) == 0 )
2220 else
2221 clear = ( cldflag == 0 )
2222 end if
2223
2224 if ( clear ) then
2225
2226 if ( iopt == 1 ) then
2227 ptop_bt(:) = min ( plev(nlev), p0 )
2228 ptop_rd(:) = min ( plev(nlev), p0 )
2229 else if ( iopt == 2 ) then
2230 ptop_bt(:) = 0.d0
2231 ptop_rd(:) = 0.d0
2232 end if
2233
2234 ntop_bt(:) = 1
2235 ntop_rd(:) = 1
2236
2237 lev_start = max ( lev_start , i10 )
2238
2239 return
2240
2241 end if
2242
2243 allocate ( ht(nlev) )
2244 nch = size( ilist)
2245
2246 channels: do jch = 1, nch
2247
2248 jc = ilist(jch)
2249
2250 ! missing channel ... yes it can happen
2251 if (jc == -1) cycle channels
2252 ! gross check failure
2253 if ( present(rejflag_opt) ) then
2254 if ( rejflag_opt(jc,9) == 1 ) cycle channels
2255 else
2256 if ( btObs(jc) < 150.d0 .or. btObs(jc) > 350.d0) cycle channels
2257 end if
2258 ! no clouds if observed radiance warmer than clear estimate
2259
2260 if ( radObs(jc) > rcal(jc) ) then
2261
2262 if ( iopt == 1 ) then
2263 ptop_bt(jc) = min ( plev(nlev), p0 )
2264 ptop_rd(jc) = min ( plev(nlev), p0 )
2265 else if ( iopt == 2 ) then
2266 ptop_bt(jc) = 0.d0
2267 ptop_rd(jc) = 0.d0
2268 end if
2269
2270 ntop_bt(jc) = 1
2271 ntop_rd(jc) = 1
2272
2273 cycle channels
2274
2275 end if
2276
2277 ! cloudy
2278
2279 if ( present(rejflag_opt)) then
2280 cloudy = ( rejflag_opt(jc,11) == 1 .and. rejflag_opt(jc,23) == 1 )
2281 else
2282 cloudy = ( cldflag == 1 )
2283 end if
2284
2285 if ( cloudy ) then
2286
2287 if ( iopt == 1 ) then
2288
2289 if ( ihgt == 0 .or. ihgt == 2 ) then
2290 call get_top ( ht,nht, btObs(jc), tt, plev, lev_start, iopt )
2291 itop = 1
2292 if ( nht >= 2 ) itop = 2
2293 ptop_bt(jc) = min ( ht(itop), p0 )
2294 ntop_bt(jc) = nht
2295 end if
2296
2297 if ( ihgt == 1 .or. ihgt == 2 ) then
2298 call get_top ( ht, nht, radObs(jc), cloudyRadiance(jc,:), plev, lev_start, iopt )
2299 itop = 1
2300 if ( nht >= 2 ) itop = 2
2301 ptop_rd(jc) = min ( ht(itop), p0 )
2302 ntop_rd(jc) = nht
2303 end if
2304
2305 else if ( iopt == 2 ) then
2306
2307 if ( ihgt == 0 .or. ihgt == 2 ) then
2308 call get_top ( ht,nht, btObs(jc),tt,height,lev_start,iopt)
2309 itop = 1
2310 if ( nht >= 2 ) itop = 2
2311 ptop_bt(jc) = max ( ht(itop), 0.d0 )
2312 ntop_bt(jc) = nht
2313 end if
2314
2315 if ( ihgt == 1 .or. ihgt == 2 ) then
2316 call get_top ( ht,nht, radObs(jc),cloudyRadiance(jc,:),height,lev_start,iopt)
2317 itop = 1
2318 if ( nht >= 2 ) itop = 2
2319 ptop_rd(jc) = max ( ht(itop), 0.d0 )
2320 ntop_rd(jc) = nht
2321 end if
2322
2323 end if
2324
2325 end if
2326
2327 end do channels
2328
2329 deallocate ( ht )
2330
2331 end subroutine cloud_top
2332
2333 !--------------------------------------------------------------------------
2334 ! get_top
2335 !--------------------------------------------------------------------------
2336 subroutine get_top (ht, nht, bt, tt ,pp, lev_start, iopt)
2337 !
2338 !:Purpose: Computation of cloud top height and number of possible heights.
2339 !
2340 implicit none
2341
2342 ! Arguments:
2343 real(8), intent(out) :: ht(:) ! Cloud top height in hpa or meters (iopt = 1 or 2)
2344 integer, intent(out) :: nht ! Number of possible cloud height solutions
2345 real(8), intent(in) :: bt ! Observed brightness temperatures (deg k) or radiance (mw/m2/sr/cm-1)
2346 real(8), intent(in) :: tt(:) ! Temperature profile (deg k) or computed cloud radiance from each level to top
2347 real(8), intent(in) :: pp(size(tt)) ! Pressure (hpa) or heights (m) profile (iopt=1 or 2)
2348 integer, intent(inout) :: lev_start ! Level to start iteration (ideally tropopause, if <= 0, search & start at coldest level)
2349 integer, intent(in) :: iopt ! Height units in hpa (1) or in meters (2)
2350
2351 ! Locals:
2352 integer :: i, im(1), i10, nlev
2353 real(8),allocatable :: logp(:)
2354 real(8) :: dt, a, b
2355
2356 ht(:) = -1.
2357
2358 im = lev_start
2359
2360 nlev = size( tt )
2361
2362 if ( lev_start <= 0 ) then
2363
2364 ! Search index im(1) where tt is minimum
2365 im = minloc ( tt )
2366
2367 i10 = -1
2368 do i = 2, nlev
2369 if (pp(i-1) <= 100.d0 .and. pp(i) > 100.d0) then
2370 I10 = I
2371 exit
2372 end if
2373 end do
2374
2375 lev_start = im(1)
2376
2377 if ( im(1) == nlev ) then
2378 lev_start = max(lev_start,i10)
2379 nht = 1
2380 ht(1) = pp(nlev)
2381 return
2382 end if
2383
2384 end if
2385
2386 if (iopt == 1) then
2387 allocate ( logp(nlev) )
2388 logp(:) = log(pp(:))
2389 end if
2390
2391 nht = 0
2392
2393 do I = im(1), nlev - 1
2394 dt = tt(I + 1) - tt(I) + 1.D-12
2395 if ( bt > tt(I) .and. bt <= tt(I + 1) ) then
2396
2397 nht = nht + 1
2398
2399 if (iopt == 1) then
2400 a = logp(I) + (logp(I + 1) - logp(I)) / dt * ( bt - tt(I))
2401 ht(nht) = exp(a)
2402 end if
2403
2404 if (iopt == 2) then
2405 b = pp(I) + (pp(I+1) - pp(I)) / dt * (bt - tt(I))
2406 ht(nht) = b
2407 end if
2408
2409 else if ( bt >= tt(I+1) .and. bt < tt(I) ) then
2410
2411 nht = nht + 1
2412
2413 if (iopt == 1) then
2414 a = logp(I + 1)- (logp(I + 1)-logp(I)) / dt * (tt(I + 1) - bt)
2415 ht(nht) = exp(A)
2416 end if
2417
2418 if(iopt == 2) then
2419 b = pp(I + 1)- (pp(I + 1) - pp(I)) / dt * (tt(I + 1) - bt)
2420 ht(nht) = b
2421 end if
2422
2423 end if
2424 end do
2425
2426
2427 if ( nht == 0 .and. bt < tt(im(1)) ) then
2428 nht = 1
2429 ht(1) = pp(im(1))
2430 else if ( nht == 0 .and. bt > tt(nlev) ) then
2431 nht = 1
2432 ht(1) = pp(nlev)
2433 end if
2434
2435 if (iopt==1) deallocate ( logp )
2436
2437 end subroutine get_top
2438
2439 !--------------------------------------------------------------------------
2440 ! get_avhrr_emiss
2441 !--------------------------------------------------------------------------
2442 subroutine get_avhrr_emiss( iasi_surfem1, freqiasi, nchaniasi, avhrr_surfem1 )
2443 !
2444 !:Purpose: choisi l'emissivite d'un canal IASI proche pour AVHRR
2445 ! a raffiner pour prendre en compte la largeur des canaux AVHRR ??
2446 !
2447 implicit none
2448
2449 ! Arguments:
2450 real(8), intent(in) :: iasi_surfem1(nchaniasi)! IASI emissivities
2451 real(8), intent(in) :: freqiasi(nchaniasi) ! IASI wavenumbers (cm-1)
2452 integer, intent(in) :: nchaniasi ! Number of IASI channels
2453 real(8), intent(out) :: avhrr_surfem1(NIR) ! AVHRR emissivities
2454
2455 ! Locals:
2456 real(8),parameter :: freqavhrr(NIR)= (/0.2687000000D+04 , 0.9272000000D+03 , 0.8377000000D+03/)
2457 integer :: indxavhrr(NIR)
2458 integer :: i, pos(1)
2459
2460 do I=1,NIR
2461 pos = minloc ( abs (freqiasi(:) - freqavhrr(I)) )
2462 indxavhrr(i) = pos(1)
2463 end do
2464
2465 do I=1,NIR
2466 avhrr_surfem1(i) = iasi_surfem1(indxavhrr(i))
2467 end do
2468
2469 end subroutine get_avhrr_emiss
2470
2471 !--------------------------------------------------------------------------
2472 ! tovs_rttov_avhrr_for_IASI
2473 !--------------------------------------------------------------------------
2474 subroutine tovs_rttov_avhrr_for_IASI (headerIndex, surfem1_avhrr, idiasi)
2475 !
2476 !:Purpose: Computation of forward radiance with rttov_direct
2477 ! (for AVHRR).
2478 ! appel de RTTOV pour le calcul des radiances AVHRR
2479 ! (non assimilees mais necessaires au background check IASI)
2480 !
2481 implicit none
2482
2483 ! Arguments:
2484 integer, intent(in) :: headerIndex ! Location of IASI observation in TOVS structures and obSpaceData
2485 real(8), intent(in) :: surfem1_avhrr(3) ! AHVRR surface emissivities
2486 integer, intent(in) :: idiasi ! iasi (in fact METOP) number
2487
2488 ! Locals:
2489 type (rttov_profile), pointer :: profiles(:)
2490 type (rttov_chanprof) :: chanprof(3)
2491 logical :: calcemis (3)
2492 integer :: list_sensor (3),errorstatus,allocStatus
2493 integer, save :: idiasi_old=-1
2494 integer :: ich
2495 integer :: ichan_avhrr (NIR)
2496 type (rttov_transmission) :: transmission
2497 type (rttov_radiance) :: radiancedata_d
2498 type (rttov_emissivity) :: emissivity(3)
2499 integer :: nchannels
2500 integer :: nlevels, iptobs(1)
2501
2502 if (idiasi_old /= idiasi) then
2503 list_sensor(1) = 10
2504 list_sensor(2) = idiasi
2505 list_sensor(3) = 5
2506 do ich=1,nir
2507 ichan_avhrr(ich)=ich
2508 end do
2509
2510 errorstatus = 0
2511
2512 if (idiasi_old > 0) then
2513 call rttov_dealloc_coefs(errorstatus, coefs_avhrr )
2514 if ( errorstatus /= 0) then
2515 write(*,*) "Probleme dans rttov_dealloc_coefs !"
2516 call utl_abort("tovs_rttov_avhrr_for_IASI")
2517 end if
2518 end if
2519
2520 call rttov_read_coefs ( errorstatus, &! out
2521 coefs_avhrr, &! out
2522 tvs_opts(1), &! in
2523 channels=ichan_avhrr, &! in
2524 instrument=list_sensor ) ! in
2525
2526 if ( errorstatus /= 0) then
2527 write(*,*) "Probleme dans rttov_read_coefs !"
2528 call utl_abort("tovs_rttov_avhrr_for_IASI")
2529 end if
2530
2531 idiasi_old = idiasi
2532
2533 end if
2534
2535 call tvs_getProfile(profiles, 'nl')
2536
2537 iptobs(1) = headerIndex
2538 nlevels = profiles(headerIndex)% nlevels
2539
2540 nchannels = NIR
2541
2542 calcemis(:) = .false.
2543 emissivity(1:3)%emis_in = surfem1_avhrr(1:3)
2544 ! Build the list of channels/profiles indices
2545
2546 do ich = 1, nchannels
2547 chanprof(ich) % prof = 1
2548 chanprof(ich) % chan = ich
2549 end do
2550
2551 allocStatus = 0
2552 call rttov_alloc_direct( &
2553 allocStatus, &
2554 asw=1, &
2555 nprofiles=1, & ! (not used)
2556 nchanprof=nchannels, &
2557 nlevels=nlevels, &
2558 opts=tvs_opts(1), &
2559 coefs=coefs_avhrr, &
2560 transmission=transmission, &
2561 radiance=radiancedata_d, &
2562 init=.true.)
2563
2564 if (allocStatus /= 0) then
2565 write(*,*) "Memory allocation error"
2566 call utl_abort('tovs_rttov_avhrr_for_IASI')
2567 end if
2568
2569
2570 call rttov_direct( &
2571 errorstatus, & ! out
2572 chanprof, & ! in
2573 tvs_opts(1), & ! in
2574 profiles(iptobs(:)), & ! in
2575 coefs_avhrr, & ! in
2576 transmission, & ! inout
2577 radiancedata_d, & ! out
2578 calcemis=calcemis, & ! in
2579 emissivity=emissivity) ! inout
2580
2581 avhrr_bgck(headerIndex)% radclearcalc(NVIS+1:NVIS+NIR) = radiancedata_d % clear(1:NIR)
2582 avhrr_bgck(headerIndex)% tbclearcalc(NVIS+1:NVIS+NIR) = radiancedata_d % bt(1:NIR)
2583 allocate( avhrr_bgck(headerIndex)% radovcalc(nlevels-1,NVIS+1:NVIS+NIR) )
2584 avhrr_bgck(headerIndex)% radovcalc(1:nlevels-1,NVIS+1:NVIS+NIR) = radiancedata_d % overcast(1:nlevels-1,1:NIR)
2585 avhrr_bgck(headerIndex)% emiss(NVIS+1:NVIS+NIR) = emissivity(1:NIR)%emis_out
2586 avhrr_bgck(headerIndex)% transmsurf(NVIS+1:NVIS+NIR) = transmission% tau_total(1:NIR)
2587
2588
2589 call rttov_alloc_direct( &
2590 allocStatus, &
2591 asw=0, &
2592 nprofiles=1, & ! (not used)
2593 nchanprof=nchannels, &
2594 nlevels=nlevels, &
2595 opts=tvs_opts(1), &
2596 coefs=coefs_avhrr, &
2597 transmission=transmission, &
2598 radiance=radiancedata_d )
2599
2600 if (allocStatus /= 0) then
2601 write(*,*) "Memory deallocation error"
2602 call utl_abort('tovs_rttov_avhrr_for_IASI')
2603 end if
2604
2605 nullify(profiles)
2606
2607 end subroutine tovs_rttov_avhrr_for_IASI
2608
2609 !--------------------------------------------------------------------------
2610 ! cor_albedo
2611 !--------------------------------------------------------------------------
2612 subroutine cor_albedo(delta, scos)
2613 !
2614 !:Purpose: ce sous-programme calcule un facteur de correction
2615 ! pour l'albedo a partir du cosinus de l'angle solaire.
2616 !
2617 implicit none
2618
2619 ! Arguments:
2620 real(8), intent(in) :: scos ! Cosine of solar zenith angle
2621 real(8), intent(out) :: delta ! Correction factor
2622
2623 ! Locals:
2624 integer i1, i2, ierr
2625 real(8) x1, x2, g1, g2, a, b
2626 real(8),parameter :: s(11)=(/00.00d0, 18.19d0, 31.79d0, 41.41d0, 49.46d0, &
2627 56.63d0, 63.26d0, 69.51d0, 75.52d0, 81.37d0, 87.13d0/)
2628
2629 i1 = 12 - ( scos + 0.05d0) * 10.d0
2630 i2 = i1 + 1
2631 i1 = min(i1,11)
2632 i2 = min(i2,11)
2633 x1 = cos ( s(I1) * MPC_RADIANS_PER_DEGREE_R8 )
2634 x2 = cos ( s(I2) * MPC_RADIANS_PER_DEGREE_R8 )
2635 g1 = drcld(i1)
2636 g2 = drcld(i2)
2637 if (i1 == i2) then
2638 delta =g1
2639 else
2640 call lineq ( x1, x2, g1, g2, a, b, ierr )
2641 delta = a * scos + b
2642 end if
2643
2644 end subroutine cor_albedo
2645
2646 !--------------------------------------------------------------------------
2647 ! drcld
2648 !--------------------------------------------------------------------------
2649 real(8) function drcld(iz)
2650 !
2651 !:Purpose: Generaliser pour toutes les plateformes satellitaires.
2652 ! Ce sous-programme calcule la normalisation due
2653 ! a l'angle zenith solaire selon
2654 ! MINNIS-HARRISSON (COURBE FIG 7), P1038,JCAM 84.
2655 !
2656 !:Output: facteur de normalisation
2657 !
2658 implicit none
2659
2660 ! Arguments:
2661 integer, intent(in) :: iz ! Index for Sun angle bin
2662
2663 ! Locals:
2664 real(8),parameter :: drf(11)=(/1.000d0, 1.002d0, 1.042d0, 1.092d0, 1.178d0, 1.286d0, &
2665 1.420d0, 1.546d0, 1.710d0, 1.870d0, 2.050d0/)
2666
2667 drcld = drf (iz)
2668
2669 end function drcld
2670
2671
2672 !--------------------------------------------------------------------------
2673 ! visocn
2674 !--------------------------------------------------------------------------
2675 subroutine visocn(sz, satz, rz, anisot, zlamb, zcloud, ierr)
2676 !
2677 !:Purpose: This routine provides the corrective factors for the anisotropy
2678 ! of reflectance over clear ocean.
2679 !
2680 !
2681 !:Notes: Obtained from dr pat minnis,langley , and based on the work
2682 ! of minnis and harrisson,jcam 1984,p993.
2683 ! the routine is a look up table along with interpolation on the
2684 ! three angles.
2685 !
2686 implicit none
2687
2688 ! Arguments:
2689 real(8), intent(in) :: sz ! sun zenith angle in degrees (0 to 90)
2690 real(8), intent(in) :: satz ! satellite zenith angle (0 to 90)
2691 real(8), intent(in) :: rz ! relative angle in degrees (0-180) with 0 as backscattering and 180 as forward scattering
2692 real(8), intent(out) :: anisot ! anisotropic corrective factor (khi in minnis-harrisson)
2693 real(8), intent(out) :: zlamb ! corrective factor for lambertian reflectance (Ocean surface)
2694 real(8), intent(out) :: zcloud ! Same as zlamb but for cloud surface
2695 integer, intent(out) :: ierr ! Error code (0=ok; -1=problem with interpolation)
2696
2697 ! Locals:
2698 integer i1, i2, j1, j2, k1, k2, l, i, n, m, j, k
2699 real(8) cc, d1, d2, slope, intercept, x1, x2
2700 real(8) g1, g2, da(2), dd(2)
2701 real(8), parameter :: s(11)=(/0.0d0,18.19d0,31.79d0,41.41d0,49.46d0,56.63d0, &
2702 63.26d0,69.51d0,75.52d0,81.37d0,87.13d0/)
2703 real(8), parameter :: r(13)=(/0.0d0, 15.0d0, 30.0d0, 45.0d0, 60.0d0, 75.0d0, 90.0d0, &
2704 105.0d0, 120.0d0, 135.0d0, 150.0d0, 165.0d0, 180.0d0/)
2705 real(8), parameter :: v(10)=(/0.0d0, 10.0d0, 20.0d0, 30.0d0, 40.0d0, 50.0d0, 60.0d0, &
2706 70.0d0, 80.0d0, 90.0d0/)
2707 real(8) vnorm(11,10,13)
2708
2709 data ((vnorm(1,j,k),j=1,10),k=1,13)/ &
2710 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2711 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2712 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2713 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2714 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2715 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2716 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2717 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2718 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2719 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2720 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2721 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0, &
2722 2.668d0,2.210d0,1.105d0,0.979d0,0.810d0,0.735d0,0.785d0,0.979d0,1.092d0,1.174d0/
2723
2724 data ((vnorm(2,j,k),j=1,10),k=1,13)/ &
2725 1.154d0, .960d0, .896d0, .818d0, .748d0, .825d0, .922d0,1.018d0,1.179d0,1.334d0, &
2726 1.154d0, .954d0, .838d0, .799d0, .735d0, .786d0, .883d0, .960d0,1.128d0,1.250d0, &
2727 1.514d0, .973d0, .825d0, .786d0, .722d0, .754d0,0.838d0,0.922d0,1.063d0,1.160d0, &
2728 1.514d0,0.967d0,0.864d0,0.818d0,0.715d0,0.728d0,0.793d0,0.876d0,1.005d0,1.102d0, &
2729 1.514d0,0.967d0,0.896d0,0.889d0,0.702d0,0.696d0,0.773d0,0.851d0,0.954d0,1.038d0, &
2730 1.514d0,1.070d0,0.986d0,0.922d0,0.677d0,0.696d0,0.754d0,0.838d0,0.922d0,1.012d0, &
2731 1.514d0,1.270d0,0.967d0,0.870d0,0.677d0,0.664d0,0.709d0,0.773d0,0.857d0,0.954d0, &
2732 1.514d0,1.495d0,1.166d0,0.960d0,0.683d0,0.690d0,0.728d0,0.806d0,0.896d0,0.999d0, &
2733 1.514d0,1.959d0,1.534d0,1.025d0,0.973d0,0.709d0,0.754d0,0.857d0,0.954d0,1.050d0, &
2734 1.514d0,2.165d0,2.165d0,1.270d0,1.038d0,0.760d0,0.812d0,0.902d0,1.012d0,1.115d0, &
2735 1.514d0,2.275d0,2.262d0,1.688d0,1.115d0,0.780d0,0.857d0,0.954d0,1.070d0,1.173d0, &
2736 1.514d0,2.326d0,2.520d0,2.172d0,1.257d0,0.812d0,0.883d0,1.005d0,1.108d0,1.212d0, &
2737 1.514d0,2.359d0,2.951d0,2.255d0,1.411d0,0.980d0,0.915d0,1.050d0,1.160d0,1.295d0/
2738
2739 data ((vnorm(3,j,k),j=1,10),k=1,13)/ &
2740 0.897d0,0.792d0,0.765d0,0.765d0,0.778d0,0.897d0,0.996d0,1.095d0,1.306d0,1.431d0, &
2741 0.897d0,0.712d0,0.739d0,0.745d0,0.765d0,0.891d0,0.970d0,1.069d0,1.214d0,1.359d0, &
2742 0.897d0,0.666d0,0.699d0,0.745d0,0.759d0,0.811d0,0.917d0,1.042d0,1.148d0,1.306d0, &
2743 0.897d0,0.646d0,0.693d0,0.739d0,0.693d0,0.752d0,0.858d0,0.989d0,1.102d0,1.234d0, &
2744 0.897d0,0.686d0,0.679d0,0.726d0,0.679d0,0.693d0,0.792d0,0.924d0,1.049d0,1.154d0, &
2745 0.897d0,0.660d0,0.673d0,0.693d0,0.646d0,0.660d0,0.759d0,0.858d0,1.003d0,1.102d0, &
2746 0.897d0,0.673d0,0.765d0,0.792d0,0.712d0,0.600d0,0.699d0,0.811d0,0.963d0,1.055d0, &
2747 0.897d0,0.706d0,0.772d0,0.917d0,0.904d0,0.613d0,0.726d0,0.858d0,1.055d0,1.121d0, &
2748 0.897d0,0.825d0,0.924d0,0.996d0,0.989d0,0.686d0,0.778d0,0.937d0,1.115d0,1.181d0, &
2749 0.897d0,1.036d0,1.253d0,1.286d0,1.260d0,0.778d0,0.858d0,0.996d0,1.181d0,1.260d0, &
2750 0.897d0,1.201d0,1.788d0,1.986d0,1.827d0,0.884d0,0.851d0,1.062d0,1.227d0,1.333d0, &
2751 0.897d0,1.530d0,2.249d0,2.546d0,2.381d0,1.352d0,0.891d0,1.108d0,1.286d0,1.405d0, &
2752 0.897d0,1.854d0,2.401d0,3.325d0,2.559d0,1.590d0,0.937d0,1.168d0,1.214d0,1.425d0/
2753
2754 data ((vnorm(4,j,k),j=1,10),k=1,13)/ &
2755 0.752d0,0.800d0,0.745d0,0.717d0,0.759d0,0.891d0,1.149d0,1.309d0,1.469d0,1.650d0, &
2756 0.752d0,0.773d0,0.717d0,0.703d0,0.752d0,0.835d0,1.065d0,1.246d0,1.406d0,1.552d0, &
2757 0.752d0,0.731d0,0.689d0,0.703d0,0.745d0,0.814d0,0.988d0,1.176d0,1.323d0,1.476d0, &
2758 0.752d0,0.689d0,0.675d0,0.654d0,0.696d0,0.752d0,0.940d0,1.100d0,1.246d0,1.378d0, &
2759 0.752d0,0.675d0,0.661d0,0.633d0,0.668d0,0.717d0,0.877d0,1.030d0,1.176d0,1.309d0, &
2760 0.752d0,0.647d0,0.640d0,0.620d0,0.613d0,0.682d0,0.814d0,0.947d0,1.107d0,1.232d0, &
2761 0.752d0,0.633d0,0.620d0,0.613d0,0.606d0,0.640d0,0.773d0,0.898d0,1.044d0,1.162d0, &
2762 0.752d0,0.626d0,0.626d0,0.626d0,0.620d0,0.654d0,0.821d0,0.947d0,1.128d0,1.225d0, &
2763 0.752d0,0.633d0,0.633d0,0.633d0,0.647d0,0.675d0,0.877d0,1.009d0,1.183d0,1.274d0, &
2764 0.752d0,0.682d0,0.717d0,0.961d0,1.023d0,0.968d0,0.940d0,1.142d0,1.274d0,1.413d0, &
2765 0.752d0,0.856d0,1.037d0,1.434d0,1.594d0,1.441d0,1.044d0,1.225d0,1.323d0,1.545d0, &
2766 0.752d0,1.044d0,1.295d0,2.207d0,1.610d0,2.311d0,1.385d0,1.274d0,1.441d0,1.636d0, &
2767 0.752d0,1.079d0,1.524d0,2.541d0,3.564d0,3.014d0,1.942d0,1.462d0,1.552d0,1.726d0/
2768
2769 data ((vnorm(5,j,k),j=1,10),k=1,13)/ &
2770 0.552d0,0.588d0,0.617d0,0.638d0,0.724d0,0.860d0,1.133d0,1.362d0,1.556d0,1.678d0, &
2771 0.552d0,0.581d0,0.602d0,0.617d0,0.652d0,0.803d0,1.075d0,1.326d0,1.484d0,1.592d0, &
2772 0.552d0,0.559d0,0.588d0,0.595d0,0.617d0,0.731d0,1.018d0,1.283d0,1.412d0,1.527d0, &
2773 0.552d0,0.531d0,0.538d0,0.574d0,0.595d0,0.710d0,0.946d0,1.240d0,1.341d0,1.463d0, &
2774 0.552d0,0.516d0,0.523d0,0.552d0,0.559d0,0.695d0,0.911d0,1.226d0,1.291d0,1.412d0, &
2775 0.552d0,0.516d0,0.523d0,0.538d0,0.538d0,0.652d0,0.882d0,1.154d0,1.240d0,1.348d0, &
2776 0.552d0,0.516d0,0.523d0,0.538d0,0.523d0,0.595d0,0.774d0,1.075d0,1.169d0,1.269d0, &
2777 0.552d0,0.531d0,0.545d0,0.552d0,0.566d0,0.609d0,0.817d0,1.140d0,1.248d0,1.369d0, &
2778 0.552d0,0.538d0,0.545d0,0.566d0,0.581d0,0.645d0,0.911d0,1.240d0,1.319d0,1.441d0, &
2779 0.552d0,0.566d0,0.552d0,0.574d0,0.710d0,0.839d0,0.982d0,1.298d0,1.391d0,2.323d0, &
2780 0.552d0,0.566d0,0.559d0,0.710d0,1.147d0,1.176d0,1.040d0,1.348d0,1.671d0,2.674d0, &
2781 0.552d0,0.588d0,1.133d0,1.355d0,2.194d0,2.803d0,2.201d0,2.459d0,2.904d0,3.126d0, &
2782 0.552d0,0.710d0,1.341d0,1.757d0,3.026d0,3.900d0,4.445d0,4.503d0,4.445d0,4.503d0/
2783
2784 data ((vnorm(6,j,k),j=1,10),k=1,13)/ &
2785 0.551d0,0.627d0,0.665d0,0.734d0,0.826d0,0.971d0,1.231d0,1.537d0,1.721d0,1.866d0, &
2786 0.551d0,0.604d0,0.619d0,0.665d0,0.765d0,0.895d0,1.185d0,1.476d0,1.568d0,1.652d0, &
2787 0.551d0,0.597d0,0.604d0,0.619d0,0.734d0,0.849d0,1.101d0,1.346d0,1.453d0,1.568d0, &
2788 0.551d0,0.581d0,0.589d0,0.597d0,0.665d0,0.795d0,1.032d0,1.262d0,1.346d0,1.445d0, &
2789 0.551d0,0.558d0,0.558d0,0.566d0,0.612d0,0.727d0,0.987d0,1.201d0,1.262d0,1.399d0, &
2790 0.551d0,0.505d0,0.505d0,0.512d0,0.566d0,0.696d0,0.925d0,1.117d0,1.185d0,1.308d0, &
2791 0.551d0,0.474d0,0.497d0,0.512d0,0.535d0,0.673d0,0.864d0,1.048d0,1.124d0,1.216d0, &
2792 0.551d0,0.497d0,0.505d0,0.520d0,0.551d0,0.681d0,0.902d0,1.124d0,1.201d0,1.323d0, &
2793 0.551d0,0.535d0,0.535d0,0.551d0,0.566d0,0.711d0,1.017d0,1.201d0,1.269d0,1.422d0, &
2794 0.551d0,0.535d0,0.543d0,0.558d0,0.704d0,1.193d0,1.247d0,1.285d0,1.346d0,1.950d0, &
2795 0.551d0,0.543d0,0.551d0,0.581d0,0.994d0,1.545d0,1.583d0,1.354d0,2.019d0,2.883d0, &
2796 0.551d0,0.566d0,0.612d0,0.788d0,1.468d0,2.233d0,2.340d0,2.531d0,2.983d0,3.365d0, &
2797 0.551d0,0.658d0,0.665d0,1.101d0,2.134d0,3.120d0,4.221d0,4.856d0,4.956d0,5.613d0/
2798
2799 data ((vnorm(7,j,k),j=1,10),k=1,13)/ &
2800 0.545d0,0.606d0,0.683d0,0.744d0,0.798d0,0.990d0,1.228d0,1.704d0,1.850d0,2.049d0, &
2801 0.545d0,0.576d0,0.583d0,0.714d0,0.783d0,0.952d0,1.144d0,1.573d0,1.758d0,1.888d0, &
2802 0.545d0,0.560d0,0.568d0,0.629d0,0.744d0,0.875d0,1.105d0,1.504d0,1.642d0,1.788d0, &
2803 0.545d0,0.553d0,0.560d0,0.599d0,0.629d0,0.791d0,1.028d0,1.420d0,1.527d0,1.696d0, &
2804 0.545d0,0.545d0,0.553d0,0.599d0,0.606d0,0.714d0,0.990d0,1.335d0,1.451d0,1.581d0, &
2805 0.545d0,0.530d0,0.537d0,0.568d0,0.583d0,0.683d0,0.890d0,1.243d0,1.351d0,1.489d0, &
2806 0.545d0,0.491d0,0.499d0,0.507d0,0.576d0,0.622d0,0.791d0,1.182d0,1.282d0,1.389d0, &
2807 0.545d0,0.507d0,0.514d0,0.507d0,0.576d0,0.675d0,0.890d0,1.197d0,1.328d0,1.451d0, &
2808 0.545d0,0.522d0,0.537d0,0.522d0,0.591d0,0.760d0,0.944d0,1.259d0,1.389d0,1.527d0, &
2809 0.545d0,0.537d0,0.545d0,0.553d0,0.614d0,0.906d0,1.028d0,1.389d0,1.504d0,2.533d0, &
2810 0.545d0,0.553d0,0.553d0,0.576d0,0.637d0,1.036d0,1.550d0,1.658d0,1.934d0,3.277d0, &
2811 0.545d0,0.560d0,0.568d0,0.606d0,1.174d0,1.781d0,2.563d0,3.170d0,3.791d0,4.966d0, &
2812 0.545d0,0.591d0,0.614d0,1.259d0,2.065d0,2.824d0,3.761d0,4.498d0,5.902d0,6.148d0/
2813
2814 data ((vnorm(8,j,k),j=1,10),k=1,13)/ &
2815 0.514d0,0.539d0,0.596d0,0.694d0,0.832d0,1.004d0,1.444d0,1.869d0,2.203d0,2.538d0, &
2816 0.514d0,0.539d0,0.571d0,0.645d0,0.751d0,0.906d0,1.387d0,1.779d0,2.056d0,2.317d0, &
2817 0.514d0,0.547d0,0.555d0,0.612d0,0.702d0,0.824d0,1.281d0,1.681d0,1.934d0,2.203d0, &
2818 0.514d0,0.539d0,0.555d0,0.588d0,0.653d0,0.743d0,1.028d0,1.404d0,1.624d0,2.024d0, &
2819 0.514d0,0.539d0,0.547d0,0.555d0,0.588d0,0.710d0,0.889d0,1.191d0,1.420d0,1.820d0, &
2820 0.514d0,0.522d0,0.522d0,0.539d0,0.563d0,0.710d0,0.849d0,1.044d0,1.208d0,1.534d0, &
2821 0.514d0,0.481d0,0.506d0,0.514d0,0.539d0,0.694d0,0.824d0,1.028d0,1.200d0,1.371d0, &
2822 0.514d0,0.481d0,0.514d0,0.547d0,0.563d0,0.702d0,0.898d0,1.134d0,1.297d0,1.501d0, &
2823 0.514d0,0.490d0,0.514d0,0.555d0,0.588d0,0.726d0,0.955d0,1.265d0,1.379d0,1.648d0, &
2824 0.514d0,0.547d0,0.547d0,0.571d0,0.604d0,0.767d0,1.036d0,1.355d0,1.550d0,3.142d0, &
2825 0.514d0,0.563d0,0.579d0,0.604d0,0.612d0,0.832d0,1.909d0,2.848d0,3.917d0,4.790d0, &
2826 0.514d0,0.522d0,0.563d0,0.677d0,0.767d0,1.420d0,2.040d0,3.158d0,4.863d0,6.291d0, &
2827 0.514d0,0.588d0,0.588d0,0.612d0,0.824d0,2.032d0,3.109d0,4.969d0,6.846d0,7.695d0/
2828
2829 data ((vnorm(9,j,k),j=1,10),k=1,13)/ &
2830 0.572d0,0.608d0,0.679d0,0.751d0,0.831d0,1.001d0,1.377d0,1.913d0,2.512d0,2.879d0, &
2831 0.572d0,0.572d0,0.608d0,0.679d0,0.760d0,0.930d0,1.243d0,1.707d0,2.369d0,2.700d0, &
2832 0.572d0,0.563d0,0.590d0,0.644d0,0.706d0,0.831d0,1.171d0,1.618d0,2.190d0,2.378d0, &
2833 0.572d0,0.554d0,0.563d0,0.599d0,0.662d0,0.760d0,1.010d0,1.502d0,2.011d0,2.235d0, &
2834 0.572d0,0.545d0,0.563d0,0.590d0,0.626d0,0.715d0,0.885d0,1.323d0,1.815d0,2.119d0, &
2835 0.572d0,0.527d0,0.554d0,0.572d0,0.608d0,0.670d0,0.724d0,1.144d0,1.618d0,1.868d0, &
2836 0.572d0,0.545d0,0.572d0,0.572d0,0.599d0,0.662d0,0.724d0,1.117d0,1.484d0,1.761d0, &
2837 0.572d0,0.554d0,0.590d0,0.599d0,0.608d0,0.679d0,0.760d0,1.216d0,1.582d0,1.922d0, &
2838 0.572d0,0.572d0,0.599d0,0.608d0,0.635d0,0.715d0,0.822d0,1.377d0,1.707d0,2.056d0, &
2839 0.572d0,0.590d0,0.608d0,0.635d0,0.662d0,0.742d0,0.912d0,1.529d0,3.075d0,4.693d0, &
2840 0.572d0,0.590d0,0.626d0,0.644d0,0.670d0,0.760d0,1.109d0,1.564d0,3.111d0,4.702d0, &
2841 0.572d0,0.599d0,0.644d0,0.662d0,0.688d0,0.822d0,1.788d0,2.816d0,5.346d0,7.295d0, &
2842 0.572d0,0.608d0,0.662d0,0.670d0,0.715d0,1.851d0,3.227d0,4.810d0,6.669d0,9.557d0/
2843
2844 data ((vnorm(10,j,k),j=1,10),k=1,13)/ &
2845 0.552d0,0.606d0,0.639d0,0.671d0,0.704d0,0.899d0,1.223d0,2.479d0,3.194d0,3.573d0, &
2846 0.552d0,0.574d0,0.606d0,0.628d0,0.682d0,0.855d0,1.148d0,2.339d0,2.642d0,3.378d0, &
2847 0.552d0,0.563d0,0.552d0,0.595d0,0.639d0,0.834d0,1.061d0,2.014d0,2.404d0,2.891d0, &
2848 0.552d0,0.563d0,0.509d0,0.552d0,0.628d0,0.801d0,0.985d0,1.689d0,2.176d0,2.653d0, &
2849 0.552d0,0.574d0,0.509d0,0.520d0,0.585d0,0.747d0,0.888d0,1.332d0,1.970d0,2.458d0, &
2850 0.552d0,0.531d0,0.509d0,0.509d0,0.531d0,0.682d0,0.801d0,1.191d0,1.819d0,2.425d0, &
2851 0.552d0,0.498d0,0.498d0,0.498d0,0.520d0,0.639d0,0.747d0,1.126d0,1.711d0,2.317d0, &
2852 0.552d0,0.498d0,0.509d0,0.509d0,0.541d0,0.671d0,0.780d0,1.278d0,1.862d0,2.598d0, &
2853 0.552d0,0.498d0,0.509d0,0.520d0,0.574d0,0.693d0,0.812d0,1.602d0,2.035d0,2.793d0, &
2854 0.552d0,0.520d0,0.520d0,0.531d0,0.595d0,0.725d0,0.844d0,1.916d0,2.588d0,3.768d0, &
2855 0.552d0,0.531d0,0.541d0,0.574d0,0.628d0,0.780d0,1.039d0,2.349d0,3.313d0,5.652d0, &
2856 0.552d0,0.574d0,0.563d0,0.606d0,0.660d0,0.812d0,1.797d0,3.010d0,5.478d0,7.492d0, &
2857 0.552d0,0.650d0,0.671d0,0.704d0,0.801d0,1.029d0,2.436d0,3.465d0,7.828d0,10.578d0/
2858
2859 data ((vnorm(11,j,k),j=1,10),k=1,13)/ &
2860 0.518d0,0.576d0,0.605d0,0.633d0,0.662d0,0.864d0,1.238d0,2.620d0,3.455d0,3.887d0, &
2861 0.518d0,0.547d0,0.576d0,0.576d0,0.633d0,0.835d0,1.123d0,2.447d0,2.821d0,3.656d0, &
2862 0.518d0,0.518d0,0.518d0,0.547d0,0.605d0,0.806d0,1.036d0,2.102d0,2.533d0,3.080d0, &
2863 0.518d0,0.518d0,0.461d0,0.518d0,0.576d0,0.777d0,0.950d0,1.727d0,2.274d0,2.821d0, &
2864 0.518d0,0.547d0,0.461d0,0.489d0,0.547d0,0.720d0,0.864d0,1.353d0,2.044d0,2.591d0, &
2865 0.518d0,0.489d0,0.461d0,0.461d0,0.489d0,0.662d0,0.777d0,1.180d0,1.871d0,2.562d0, &
2866 0.518d0,0.461d0,0.461d0,0.461d0,0.489d0,0.605d0,0.720d0,1.123d0,1.756d0,2.418d0, &
2867 0.518d0,0.461d0,0.461d0,0.461d0,0.518d0,0.633d0,0.749d0,1.296d0,1.929d0,2.764d0, &
2868 0.518d0,0.461d0,0.461d0,0.489d0,0.547d0,0.662d0,0.777d0,1.641d0,2.130d0,2.994d0, &
2869 0.518d0,0.489d0,0.489d0,0.489d0,0.547d0,0.691d0,0.806d0,1.986d0,2.735d0,4.117d0, &
2870 0.518d0,0.489d0,0.489d0,0.547d0,0.576d0,0.749d0,1.008d0,2.476d0,3.599d0,6.334d0, &
2871 0.518d0,0.547d0,0.518d0,0.576d0,0.633d0,0.777d0,1.842d0,3.224d0,6.132d0,8.550d0, &
2872 0.518d0,0.605d0,0.633d0,0.662d0,0.777d0,1.008d0,2.562d0,3.771d0,8.953d0,12.293d0/
2873
2874 ! compute sun zenith bin
2875 cc = cos( sz * MPC_RADIANS_PER_DEGREE_R8)
2876 i1 = 12.d0 - (cc + 0.05d0) * 10.d0
2877 i2 = i1 + 1
2878 if (i1 >= 11) i1 = 11
2879 if (i1 == 11) i2 = i1
2880
2881 ! compute sat zenith bin
2882 j1 = int(satz / 10.d0) + 1
2883 j2 = j1 + 1
2884 if (j1 == 10) j2 = j1
2885
2886 ! compute relative azimuth bin
2887 k1 = RZ / 15.d0 + 1.d0
2888 k2 = k1 + 1
2889 if (k1 == 13) k2 = k1
2890
2891 ! interpolate
2892 ierr = 0
2893 do l=i1,i2
2894 i = l -i1 + 1
2895
2896 ! between r's for constant s
2897 do n=k1,k2
2898
2899 ! between v's for constant r and s
2900 m = n - k1 + 1
2901 d1 = vnorm(l,j1,n)
2902 d2 = vnorm(l,j2,n)
2903 if (d1 == d2) then
2904 da(m) = d1
2905 else
2906 call lineq(V(j1),V(j2),d1,d2,slope,intercept,ierr)
2907 da(m) = slope * satz + intercept
2908 end if
2909 end do
2910 if(k1 == k2) then
2911 dd(i) = da(1)
2912 else
2913 call lineq(R(k1),R(k2),da(1),da(2),slope,intercept,ierr)
2914 dd(i) = slope * RZ + intercept
2915 end if
2916 end do
2917
2918 ! between s's using result of other interpolations
2919 if(i1 == i2) then
2920 zlamb = drm(i1)
2921 zcloud = drcld(i1)
2922 anisot = dd(1)
2923 else
2924 x1 = cos(s(i1) * MPC_RADIANS_PER_DEGREE_R8)
2925 x2 = cos(s(i2) * MPC_RADIANS_PER_DEGREE_R8)
2926 call lineq(x1,x2,dd(1),dd(2),slope,intercept,ierr)
2927 anisot = slope * cc + intercept
2928 g1 = drm(i1)
2929 g2 = drm(i2)
2930 call lineq(x1,x2,g1,g2,slope,intercept,ierr)
2931 zlamb = slope * cc + intercept
2932 g1 = drcld(i1)
2933 g2 = drcld(i2)
2934 call lineq(x1,x2,g1,g2,slope,intercept,ierr)
2935 zcloud = slope * cc + intercept
2936 end if
2937
2938 if (anisot < 0.) then
2939 ierr = -1
2940 anisot = 1.d0
2941 zlamb = drm(i1)
2942 zcloud = drcld(i1)
2943 end if
2944
2945 end subroutine visocn
2946
2947 !--------------------------------------------------------------------------
2948 ! lineq
2949 !--------------------------------------------------------------------------
2950 subroutine lineq(x1, x2, y1, y2, a, b, ierr)
2951 !
2952 !:Purpose: calculate slope and intercept of a line.
2953 !
2954 implicit none
2955
2956 ! Arguments:
2957 real(8), intent(in) :: x1 ! coordinate x of point 1
2958 real(8), intent(in) :: x2 ! coordinate x of point 2
2959 real(8), intent(in) :: y1 ! coordinate y of point 1
2960 real(8), intent(in) :: y2 ! coordinate y of point 2
2961 real(8), intent(out) :: a ! slope
2962 real(8), intent(out) :: b ! intercept
2963 integer, intent(out) :: ierr ! error code (0=ok)
2964
2965 ierr = 0
2966
2967 if ( (x2 - x1) == 0.d0) then
2968 ierr = -1
2969 return
2970 end if
2971
2972 a = ( y2 - y1) / (x2 - x1)
2973 b = y1 - a * x1
2974
2975 end subroutine lineq
2976
2977
2978 !--------------------------------------------------------------------------
2979 ! drm
2980 !--------------------------------------------------------------------------
2981 real(8) function drm(iz)
2982 !
2983 !:Purpose: Normalization for sun zenith angle (lambertian)
2984 ! for ocean.
2985 !
2986 !:Outputs: normalization factor
2987 !
2988 implicit none
2989
2990 ! Arguments:
2991 integer, intent(in) :: iz ! index
2992
2993 ! Locals:
2994 real(8),parameter :: drf(11)=(/1.d0,1.0255d0,1.1197d0,1.2026d0,1.3472d0, &
2995 1.4926d0,1.8180d0,2.1980d0, 2.8180d0,3.8615d0,4.3555d0/)
2996
2997 drm = drf(IZ)
2998
2999 end function drm
3000
3001
3002end module multiIRbgck_mod
