1
2module tovsNL_mod
3 ! MODULE tovsNL_mod (prefix='tvs' category='5. Observation operators')
4 !
5 !:Purpose: Derived types, public variables and procedures related to the
6 ! nonlinear version of RTTOV
7 !
8 use rttovInterfaces_mod
9 use rttov_types, only : &
10 rttov_coefs ,&
11 rttov_fast_coef ,&
12 rttov_scatt_coef ,&
13 rttov_options ,&
14 rttov_options_scatt ,&
15 rttov_profile ,&
16 rttov_profile_cloud ,&
17 rttov_radiance ,&
18 rttov_transmission ,&
19 rttov_chanprof ,&
20 rttov_emissivity
21 use rttov_const, only : &
22 platform_name ,&
23 nplatforms ,&
24 inst_name ,&
25 ninst ,&
26 inst_id_goesim ,&
27 inst_id_gmsim ,&
28 inst_id_mtsatim ,&
29 inst_id_amsua ,&
30 inst_id_mhs ,&
31 sensor_id_mw ,&
32 sensor_id_po ,&
33 platform_id_jpss ,&
34 platform_id_himawari ,&
35 platform_id_eos ,&
36 errorstatus_success ,&
37 mair, mh2o, mo3 ,&
38 surftype_land ,&
39 surftype_seaice ,&
40 surftype_sea ,&
41 watertype_ocean_water ,&
42 ngases_max ,&
43 gas_id_mixed ,&
44 gas_unit_specconc ,&
45 interp_rochon_loglinear_wfn ,&
46 zenmax ,&
47 zenmaxv9 ,&
48 o3min ,&
49 o3max ,&
50 tmin ,&
51 tmax ,&
52 qmin ,&
53 qmax ,&
54 elevmax ,&
55 wmax ,&
56 pmin ,&
57 pmax
58 use parkind1, only : jpim, jplm
59 use rttov_fast_coef_utils_mod, only: set_pointers, set_fastcoef_level_bounds
60 use rttov_solar_refl_mod, only : rttov_refl_water_interp
61 use midasMpi_mod
62 use codtyp_mod
63 use mpi
64 use utilities_mod
65 use obsSpaceData_mod
66 use earthConstants_mod
67 use MathPhysConstants_mod
68 use ozoneClim_mod
69 use columnData_mod
70 use mod_rttov_emis_atlas
71 use verticalCoord_mod
72 use codePrecision_mod
73 use humidityLimits_mod
74 use interpolation_mod
75
76 implicit none
77 save
78 private
79
80 type surface_params
81 real(8) :: albedo ! surface albedo (0-1)
82 real(8) :: ice ! ice cover (0-1)
83 real(8) :: snow ! snow cover (0-1)
84 real(8) :: pcnt_wat ! water percentage in pixel containing profile (0-1)
85 real(8) :: pcnt_reg ! water percentage in an area around profile (0-1)
86 integer :: ltype ! surface type (1,...,20)
87 end type surface_params
88
89 ! public variables (parameters)
90 public :: tvs_maxChannelNumber, tvs_maxNumberOfChannels, tvs_maxNumberOfSensors, tvs_defaultEmissivity
91 ! public variables (non-parameters)
92 public :: tvs_nchan, tvs_ichan, tvs_lsensor, tvs_headerIndex, tvs_tovsIndex, tvs_nobtov
93 public :: tvs_nchanMpiGlobal, tvs_ichanMpiGlobal
94 public :: tvs_isReallyPresent,tvs_listSensors
95 public :: tvs_isReallyPresentMpiGlobal
96 public :: tvs_nsensors, tvs_platforms, tvs_satellites, tvs_instruments, tvs_channelOffset
97 public :: tvs_debug, tvs_satelliteName, tvs_instrumentName, tvs_useO3Climatology
98 public :: tvs_coefs, tvs_opts, tvs_transmission,tvs_emissivity
99 public :: tvs_coef_scatt, tvs_opts_scatt
100 public :: tvs_radiance, tvs_surfaceParameters
101 public :: tvs_numMWInstrumUsingCLW, tvs_numMWInstrumUsingHydrometeors
102 public :: tvs_mwInstrumUsingCLW_tl, tvs_mwInstrumUsingHydrometeors_tl
103 public :: tvs_mwAllskyAssim
104 ! public procedures
105 public :: tvs_fillProfiles, tvs_rttov, tvs_printDetailledOmfStatistics, tvs_allocTransmission, tvs_cleanup
106 public :: tvs_deallocateProfilesNlTlAd
107 public :: tvs_setupAlloc,tvs_setup, tvs_isIdBurpTovs, tvs_isIdBurpHyperSpectral, tvs_isIdBurpInst, tvs_getAllIdBurpTovs
108 public :: tvs_isInstrumGeostationary, tvs_isNameHyperSpectral
109 public :: tvs_isNameGeostationary
110 public :: tvs_getInstrumentId, tvs_getPlatformId, tvs_mapSat, tvs_mapInstrum
111 public :: tvs_isInstrumHyperSpectral, tvs_getChanprof, tvs_countRadiances
112 public :: tvs_ChangedStypValue
113 public :: tvs_getHIREmissivities, tvs_getOtherEmissivities, tvs_rttov_read_coefs
114 public :: tvs_getLocalChannelIndexFromChannelNumber
115 public :: tvs_getMWemissivityFromAtlas, tvs_getProfile
116 public :: tvs_getCorrectedSatelliteAzimuth
117 public :: tvs_isInstrumUsingCLW, tvs_isInstrumUsingHydrometeors, tvs_getChannelNumIndexFromPPP
118 public :: tvs_isInstrumAllskyTtAssim, tvs_isInstrumAllskyHuAssim
119 ! Module parameters
120 ! units conversion from mixing ratio to ppmv and vice versa
121 real(8), parameter :: qMixratio2ppmv = (1000000.0d0 * mair) / mh2o
122 real(8), parameter :: qppmv2Mixratio = mh2o / (1000000.0d0 * mair)
123 real(8), parameter :: o3Mixratio2ppmv = (1000000.0d0 * mair) / mo3
124 real(8), parameter :: o3ppmv2Mixratio = mo3 / (1000000.0d0 * mair)
125 real(pre_obsReal), parameter :: tvs_defaultEmissivity = 0.95
126
127 integer, parameter :: tvs_maxChannelNumber = 8461 ! Max. value for channel number
128 integer, parameter :: tvs_maxNumberOfChannels = 2211 ! Max. no. of channels (for one profile/spectra)
129 integer, parameter :: tvs_maxNumberOfSensors = 100 ! Max no sensors to be used
130 integer, parameter :: tvs_nlevels = 101 ! Maximum No. of RTTOV pressure levels including 'rttov top' at 0.005 hPa
131
132 ! Module variables
133 integer, allocatable :: tvs_bodyIndexFromBtIndex(:) ! Provides the bodyIndex in ObsSpaceData based on btIndex
134 integer, allocatable :: tvs_nchan(:) ! Number of channels per instrument (local)
135 integer, allocatable :: tvs_ichan(:,:) ! List of channels per instrument (local)
136 integer, allocatable :: tvs_nchanMpiGlobal(:) ! Number of channels per instrument (global)
137 integer, allocatable :: tvs_ichanMpiGlobal(:,:) ! List of channels per instrument (global)
138 integer, allocatable :: tvs_lsensor(:) ! Sensor number for each profile
139 integer, allocatable :: tvs_headerIndex (:) ! Observation position in obsSpaceData header for each profile
140 integer, allocatable :: tvs_tovsIndex (:) ! Index in TOVS structures for each observation header in obsSpaceData
141 logical, allocatable :: tvs_isReallyPresent(:) ! Logical flag to identify instruments really assimilated (local)
142 logical, allocatable :: tvs_isReallyPresentMpiGLobal(:) ! Logical flag to identify instruments really assimilated (global)
143 integer, allocatable :: tvs_listSensors(:,:) ! Sensor list
144 type (rttov_emis_atlas_data), allocatable :: tvs_atlas(:) ! Emissivity atlases
145 type(surface_params), allocatable :: tvs_surfaceParameters(:) ! surface parameters
146 integer tvs_nobtov ! Number of tovs observations
147 integer tvs_nsensors ! Number of individual sensors.
148 integer tvs_platforms(tvs_maxNumberOfSensors) ! RTTOV platform ID's (e.g., 1=NOAA; 2=DMSP; ...)
149 integer tvs_satellites(tvs_maxNumberOfSensors) ! RTTOV satellite ID's (e.g., 1 to 16 for NOAA; ...)
150 integer tvs_instruments(tvs_maxNumberOfSensors) ! RTTOVinstrument ID's (e.g., 3=AMSU-A; 4=AMSU-B; 6=SSMI; ...)
151 integer tvs_channelOffset(tvs_maxNumberOfSensors)! BURP to RTTOV channel mapping offset
152 integer instrumentIdsUsingCLW(tvs_maxNumberOfSensors)
153 integer instrumentIdsUsingHydrometeors(tvs_maxNumberOfSensors)
154 integer tvs_numMWInstrumUsingCLW
155 integer tvs_numMWInstrumUsingHydrometeors
156 logical tvs_mwInstrumUsingCLW_tl
157 logical tvs_mwInstrumUsingHydrometeors_tl
158 logical tvs_mwAllskyAssim
159 logical :: tvs_mpiTask0ReadCoeffs
160 real(8) :: tvs_cloudScaleFactor
161 logical tvs_debug ! Logical key controlling statements to be executed while debugging TOVS only
162 logical tvs_useO3Climatology ! Determine if ozone model field or climatology is used
163 ! If ozone model field is specified, related increments will be generated in assimilation
164 logical tvs_regLimitExtrap ! use RTTOV reg_limit_extrap option
165 logical tvs_doAzimuthCorrection(tvs_maxNumberOfSensors)
166 logical tvs_isAzimuthValid(tvs_maxNumberOfSensors)
167 logical tvs_userDefinedDoAzimuthCorrection
168 logical tvs_userDefinedIsAzimuthValid
169
170 character(len=15) tvs_satelliteName(tvs_maxNumberOfSensors)
171 character(len=15) tvs_instrumentName(tvs_maxNumberOfSensors)
172 character(len=8) radiativeTransferCode ! RadiativeTransferCode : TOVS radiation model used
173 real(8), allocatable :: tvs_emissivity(:,:) ! Surface emissivities organized by profiles and channels
174 integer, parameter :: kslon=2160, kslat=1080 ! CERES file dimension in grid points
175
176 ! High resolution surface fields
177 integer :: surfaceType(kslon,kslat)
178 real(8) :: waterFraction(kslon,kslat)
179
180 ! Derived typeso
181 type(rttov_coefs), allocatable :: tvs_coefs(:) ! rttov coefficients
182 type(rttov_options), allocatable :: tvs_opts(:) ! rttov options
183 type(rttov_scatt_coef),allocatable :: tvs_coef_scatt(:) ! rttovscatt coefficients
184 type(rttov_options_scatt), allocatable :: tvs_opts_scatt(:) ! rttovscatt options
185 type(rttov_profile), target, allocatable :: tvs_profiles_nl(:) ! all profiles on trial vertical coordinate for nl obs operator
186 type(rttov_profile), target, allocatable :: tvs_profiles_tlad(:) ! all profiles on increments vertical coordinates for linearized obs. operator
187 type(rttov_radiance), allocatable :: tvs_radiance(:) ! radiances organized by profile
188 type(rttov_transmission), allocatable :: tvs_transmission(:) ! transmittances all profiles for HIR quality control
189 type(rttov_profile_cloud), target, allocatable :: tvs_cld_profiles_nl(:)! rttov scatt cloud profiles on trial vertical coordinate
190 type(rttov_profile_cloud), target, allocatable :: tvs_cld_profiles_tlad(:) ! rttov scatt cloud profiles on increment vertical coordinates
191
192 ! Namelist variables:
193 logical useUofWIREmiss ! Flag to activate use of RTTOV U of W IR emissivity Atlases
194 logical useMWEmissivityAtlas ! Flag to activate use of RTTOV built-in MW emissivity Atlases
195 integer mWAtlasId ! MW Atlas Id used when useMWEmissivityAtlas == .true. ; 1 TELSEM2, 2 CNRM atlas
196
197 integer, external :: get_max_rss
198
199contains
200
201 !--------------------------------------------------------------------------
202 ! validateRttovVariable
203 !--------------------------------------------------------------------------
204 subroutine validateRttovVariable(value, variableName, obsSpaceData, headerIndex, minimum_opt, maximum_opt)
205 !
206 !:Purpose: check variable for validity for RTTOV-13,
207 ! if invalid replace by acceptable value and reject
208 !
209 implicit none
210
211 ! Arguments:
212 real(8), intent(inout) :: value ! Value of the RTTOV input variable to check for validity
213 character(len=*), intent(in) :: variableName ! Name of the checked variable for output in listings
214 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData object
215 integer, intent(in) :: headerIndex ! Index of the observation in obsSpaceData header table
216 real(8), optional, intent(in) :: minimum_opt ! Minimum acceptable value
217 real(8), optional, intent(in) :: maximum_opt ! Maximum acceptable value
218
219 ! Locals:
220 real(8) :: minimum, maximum
221
222 if (present(minimum_opt)) then
223 minimum = minimum_opt
224 else
225 minimum = - huge(0.d0)
226 end if
227
228 if (present(maximum_opt)) then
229 maximum = maximum_opt
230 else
231 maximum = huge(0.d0)
232 end if
233
234 if ( value < minimum .or. value > maximum ) then
235 write(*,*) 'validateRttovVariable: !!! WARNING !!!'
236 write(*,*) 'validateRttovVariable: INVALID ' // trim(variableName)
237 write(*,*) 'validateRttovVariable: headerIndex ', headerIndex, " !"
238 write(*,*) 'validateRttovVariable: ', value, ' should be between ', minimum, ' and ', maximum
239 value = max(minimum, min(value, maximum) )
240 write(*,*) 'validateRttovVariable: replaced with ', value, ' !'
241 call rejectObs(obsSpaceData, headerIndex)
242 end if
243
244 end subroutine validateRttovVariable
245
246 !--------------------------------------------------------------------------
247 ! validateRttovProfile
248 !--------------------------------------------------------------------------
249 subroutine validateRttovProfile(profile, profileName, minimum, maximum, obsSpaceData, headerIndex)
250 !
251 !:Purpose: check profile for validity for RTTOV-13,
252 ! if invalid replace by acceptable value(s) and reject
253 !
254 implicit none
255
256 ! Arguments:
257 real(8), intent(inout) :: profile(:) ! Vertical profile of input variables to check for validity
258 character(len=*), intent(in) :: profileName ! Name of the checked profile for output in listings
259 real(8), intent(in) :: minimum ! Minimum acceptable value
260 real(8), intent(in) :: maximum ! Maximum acceptable value
261 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData object
262 integer, intent(in) :: headerIndex ! Index of the observation in obsSpaceData header table
263
264 ! Locals:
265 logical :: ltest(size(profile))
266 integer :: levelIndex
267
268 ltest(:) = (profile(:) > maximum .or. profile(:) < minimum)
269
270 if ( any(ltest) ) then
271 write(*,*) 'validateRttovProfile: !!! WARNING !!!'
272 write(*,*) 'validateRttovProfile: some INVALID ' // trim(profileName)
273 write(*,*) 'validateRttovProfile: headerIndex ', headerIndex, " !"
274 do levelIndex = 1, size(profile)
275 if ( ltest(levelIndex) ) then
276 write(*,*) 'validateRttovProfile: ', profile(levelIndex), &
277 'should be between ', minimum, ' and ', maximum
278 profile(levelIndex) = max(minimum, min(profile(levelIndex), maximum) )
279 write(*,*) 'validateRttovProfile: replaced with ', profile(levelIndex)
280 end if
281 end do
282 call rejectObs(obsSpaceData, headerIndex)
283 end if
284
285 end subroutine validateRttovProfile
286
287 !--------------------------------------------------------------------------
288 ! rejectObs
289 !--------------------------------------------------------------------------
290 subroutine rejectObs(obsSpaceData, headerIndex)
291 !
292 !:Purpose: reject all data corresponding to headerIndex
293 !
294 implicit none
295
296 ! Arguments:
297 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData object
298 integer, intent(in) :: headerIndex ! Index of the observation in obsSpaceData header table
299
300 ! Locals:
301 integer :: bodyIndex
302
303 call obs_set_current_body_list(obsSpaceData, headerIndex)
304 BODY:do
305 bodyIndex = obs_getBodyIndex(obsSpaceData)
306 if (bodyIndex < 0) exit BODY
307 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
308 end do BODY
309 end subroutine rejectObs
310
311 !--------------------------------------------------------------------------
312 ! tvs_setupAlloc
313 !--------------------------------------------------------------------------
314 subroutine tvs_setupAlloc(obsSpaceData)
315 !
316 !:Purpose: Memory allocation for the non linear radiative transfer model variables.
317 !
318 implicit none
319
320 ! Arguments:
321 type(struct_obs), intent(inout) :: obsSpaceData
322
323 ! Locals:
324 integer :: allocStatus(9)
325 integer :: satelliteCode, instrumentCode, iplatform, isat, instrum
326 integer :: tovsIndex, idatyp, sensorIndex
327 integer :: channelNumber, nosensor, channelIndex
328 integer :: errorStatus
329 integer :: headerIndex, bodyIndex, taskIndex
330 logical :: runObsOperatorWithClw
331 logical :: runObsOperatorWithHydrometeors
332 logical, allocatable :: logicalBuffer(:)
333 character(len=32) :: hydroTableFilename
334
335 if (tvs_nsensors == 0) return
336
337 ! 1. Determine the number of radiances to be assimilated.
338 ! Construct a list of channels for each sensor.
339 ! Construct a list of sensor number for each profile
340
341 write(*,*) 'tvs_setupAlloc: Starting'
342
343 allocStatus = 0
344 allocate (tvs_nchan(tvs_nsensors), stat= allocStatus(1))
345 allocate (tvs_ichan(tvs_maxNumberOfChannels,tvs_nsensors), stat= allocStatus(2))
346 allocate (tvs_lsensor(obs_numheader(obsSpaceData)), stat= allocStatus(3))
347 allocate (tvs_headerIndex (obs_numheader(obsSpaceData)), stat= allocStatus(4))
348 allocate (tvs_tovsIndex(obs_numheader(obsSpaceData)), stat= allocStatus(5))
349 allocate (tvs_isReallyPresent(tvs_nsensors), stat= allocStatus(6))
350 allocate (tvs_nchanMpiGlobal(tvs_nsensors), stat= allocStatus(7))
351 allocate (tvs_ichanMpiGlobal(tvs_maxNumberOfChannels,tvs_nsensors), stat= allocStatus(8))
352 allocate (tvs_isReallyPresentMpiGlobal(tvs_nsensors), stat= allocStatus(9))
353
354 call utl_checkAllocationStatus(allocStatus, ' tvs_setupAlloc')
355
356 tvs_nchan(:) = 0
357 tvs_ichan(:,:) = 0
358 tvs_isReallyPresent(:) = .true.
359 tvs_lsensor(:) = -1
360 tvs_headerIndex(:) = -1
361 tvs_tovsIndex (:) = -1
362
363 tvs_nobtov = 0
364
365 ! Loop over all header indices of the 'TO' family
366 ! Set the header list & start at the beginning of the list
367 call obs_set_current_header_list(obsSpaceData,'TO')
368 HEADER: do
369 headerIndex = obs_getHeaderIndex(obsSpaceData)
370 if (headerIndex < 0) exit HEADER
371
372 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
373
374 if ( .not. tvs_isIdBurpTovs(idatyp) ) then
375 write(*,*) 'tvs_setupAlloc: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
376 call rejectObs(obsSpaceData, headerIndex)
377 cycle HEADER ! Proceed to the next headerIndex
378 end if
379 tvs_nobtov = tvs_nobtov + 1
380
381 ! Construct list of channels for each sensor:
382 ! map burp satellite info to RTTOV platform and satellite.
383 satelliteCode = obs_headElem_i(obsSpaceData,OBS_SAT,headerIndex)
384 call tvs_mapSat(satelliteCode,iplatform,isat)
385 if (iplatform == -1) then
386 write(*,*) 'Unknown OBS_SAT !', satelliteCode
387 call utl_abort('tvs_setupAlloc')
388 end if
389 ! map burp instrument info to RTTOV instrument.
390 instrumentCode = obs_headElem_i(obsSpaceData,OBS_INS,headerIndex)
391 call tvs_mapInstrum(instrumentCode,instrum)
392 if (instrum == -1) then
393 write(*,*) 'Unknown OBS_INS !', instrumentCode
394 call utl_abort('tvs_setupAlloc')
395 end if
396 ! find sensor number for this obs.
397 nosensor =0
398 do sensorIndex = 1, tvs_nsensors
399 if ( iplatform == tvs_platforms (sensorIndex) .and. &
400 isat == tvs_satellites (sensorIndex) .and. &
401 instrum == tvs_instruments(sensorIndex) ) then
402 nosensor = sensorIndex
403 exit
404 end if
405 end do
406
407 if (nosensor > 0) then
408 tvs_lsensor(tvs_nobtov) = nosensor
409 tvs_headerIndex (tvs_nobtov) = headerIndex
410 tvs_tovsIndex (headerIndex) = tvs_nobtov
411 else
412 write(*,*) ' tvs_setupAlloc: Warning Invalid Sensor ', iplatform, isat, instrum, ' skipping ...'
413 end if
414
415 ! Loop over all body indices (still in the 'TO' family)
416 ! Set the body list & start at the beginning of the list
417 call obs_set_current_body_list(obsSpaceData, headerIndex)
418 BODY: do
419 bodyIndex = obs_getBodyIndex(obsSpaceData)
420 if (bodyIndex < 0) exit BODY
421 if (nosensor > 0) then
422 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
423 call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
424 channelNumber, channelIndex )
425 if ( channelIndex == 0 ) then
426 tvs_nchan(nosensor) = tvs_nchan(nosensor) + 1
427 tvs_ichan(tvs_nchan(nosensor),nosensor) = channelNumber
428 end if
429
430 if ( tvs_debug .and. mmpi_myid == 0 .and. &
431 trim(tvs_instrumentName(nosensor)) == 'AMSUA' ) then
432 write(*,*) 'test channelNumber:', headerIndex, bodyIndex, nosensor, &
433 tvs_satelliteName(nosensor), channelNumber, channelIndex
434 end if
435 end if
436 else
437 ! set to notAssimilated if instrument not in NAMTOV namelist
438 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
439 end if
440 end do BODY
441 end do HEADER
442
443 if ( .not. tvs_userDefinedDoAzimuthCorrection) then
444 ! tvs_doAzimuthCorrection user defined values will be overwriten by the old default values
445 do sensorIndex = 1, tvs_nsensors
446 tvs_doAzimuthCorrection(sensorIndex) = ( tvs_platforms(sensorIndex) /= platform_id_eos .and. &
447 ( tvs_instruments(sensorIndex) == inst_id_amsua .or. tvs_instruments(sensorIndex) == inst_id_mhs ) )
448 end do
449 if ( mmpi_myId == 0 ) write(*,*) ' tvs_setupAlloc: Warning tvs_doAzimuthCorrection user defined values overwriten by the old default values'
450 end if
451
452 if ( .not. tvs_userDefinedIsAzimuthValid ) then
453 ! tvs_isAzimuthValid user defined values will be overwriten by the current default values
454 do sensorIndex = 1, tvs_nsensors
455 tvs_isAzimuthValid(sensorIndex) = .not. ( tvs_isInstrumGeostationary(tvs_instruments(sensorIndex)) )
456 end do
457 if ( mmpi_myId == 0 ) write(*,*) ' tvs_setupAlloc: Warning tvs_isAzimuthValid user defined values overwriten by the current default values'
458 end if
459
460 if ( mmpi_myId == 0 ) then
461 write(*,*) ' tvs_setupAlloc: platform satellite id tvs_doAzimuthCorrection tvs_isAzimuthValid'
462 do sensorIndex = 1, tvs_nsensors
463 write(*,'(18x,a,1x,a,1x,i2,1x,L1,10x,L1)') inst_name(tvs_instruments(sensorIndex)), &
464 platform_name(tvs_platforms(sensorIndex)), tvs_satellites(sensorIndex), &
465 tvs_doAzimuthCorrection(sensorIndex), tvs_isAzimuthValid(sensorIndex)
466 end do
467 end if
468
469 ! Sort list of channels in ascending order.Also force at least one channel, if none are found.
470 do sensorIndex = 1, tvs_nsensors
471 call isort(tvs_ichan(:,sensorIndex),tvs_nchan(sensorIndex))
472 if ( tvs_nchan(sensorIndex) == 0 ) then
473 tvs_isReallyPresent ( sensorIndex ) =.false.
474 tvs_nchan(sensorIndex) = 1
475 tvs_ichan(1,sensorIndex) = 1
476 end if
477 end do
478
479 write(*,*) ' tvs_setupAlloc: tvs_nobtov = ', tvs_nobtov
480
481 do sensorIndex = 1, tvs_nsensors
482 call tvs_getCommonChannelSet(tvs_ichan(:,sensorIndex),tvs_nchanMpiGlobal(sensorIndex), tvs_ichanMpiGlobal(:,sensorIndex))
483 print *, 'sensorIndex,tvs_nchan(sensorIndex),tvs_nchanMpiGlobal(sensorIndex)', sensorIndex, tvs_nchan(sensorIndex),tvs_nchanMpiGlobal(sensorIndex)
484 end do
485
486 if (mmpi_myid ==0) then
487 allocate(logicalBuffer(mmpi_nprocs))
488 else
489 allocate(logicalBuffer(1))
490 end if
491
492 do sensorIndex = 1, tvs_nsensors
493 call RPN_COMM_gather( tvs_isReallyPresent ( sensorIndex ) , 1, 'MPI_LOGICAL', logicalBuffer, 1,'MPI_LOGICAL', 0, 'GRID', errorStatus )
494 if (mmpi_myid ==0) then
495 tvs_isReallyPresentMpiGlobal ( sensorIndex ) =.false.
496 do taskIndex=1, mmpi_nprocs
497 tvs_isReallyPresentMpiGlobal ( sensorIndex ) = tvs_isReallyPresentMpiGlobal ( sensorIndex ) .or. logicalBuffer(taskIndex)
498 end do
499 end if
500 call rpn_comm_bcast(tvs_isReallyPresentMpiGlobal ( sensorIndex ), 1, 'MPI_LOGICAL', 0, 'GRID', errorStatus )
501 end do
502
503 deallocate(logicalBuffer)
504
505 if ( tvs_debug .and. mmpi_myid == 0 ) then
506 do sensorIndex = 1, tvs_nsensors
507 write(*,*) 'sensorIndex, tvs_instrumentName(sensorIndex), tvs_satelliteName(sensorIndex)'
508 write(*,*) sensorIndex, tvs_instrumentName(sensorIndex), tvs_satelliteName(sensorIndex)
509 write(*,*) 'tvs_channelOffset(sensorIndex), tvs_nchan(sensorIndex)'
510 write(*,*) tvs_channelOffset(sensorIndex), tvs_nchan(sensorIndex)
511 write(*,*) 'tvs_ichan(1:tvs_nchan(sensorIndex),sensorIndex)'
512 write(*,*) tvs_ichan(1:tvs_nchan(sensorIndex),sensorIndex)
513 write(*,*)
514 end do
515 end if
516
517 ! 3. Initialize TOVS radiance transfer model
518
519 if ( radiativeTransferCode == 'RTTOV' ) then
520
521 write(*,'(//,10x,A)') '-rttov_setup: initializing the TOVS radiative transfer model'
522
523 allocate (tvs_coefs(tvs_nsensors) ,stat= allocStatus(1))
524 allocate (tvs_listSensors (3,tvs_nsensors) ,stat= allocStatus(2))
525 allocate (tvs_opts (tvs_nsensors) ,stat= allocStatus(3))
526 if (tvs_numMWInstrumUsingHydrometeors > 0) then
527 allocate (tvs_opts_scatt (tvs_nsensors) ,stat= allocStatus(4))
528 allocate (tvs_coef_scatt (tvs_nsensors) ,stat= allocStatus(5))
529 end if
530 call utl_checkAllocationStatus(allocStatus(1:5), ' tvs_setupAlloc before rttov initialization')
531
532 do sensorIndex=1, tvs_nsensors
533 tvs_listSensors(1,sensorIndex) = tvs_platforms (sensorIndex)
534 tvs_listSensors(2,sensorIndex) = tvs_satellites (sensorIndex)
535 tvs_listSensors(3,sensorIndex) = tvs_instruments(sensorIndex)
536
537 runObsOperatorWithClw = (tvs_numMWInstrumUsingCLW /= 0 .and. &
538 tvs_isInstrumUsingCLW(tvs_instruments(sensorIndex)))
539 runObsOperatorWithHydrometeors = (tvs_numMWInstrumUsingHydrometeors /= 0 .and. &
540 tvs_isInstrumUsingHydrometeors(tvs_instruments(sensorIndex)))
541
542 !< General configuration options
543 tvs_opts(sensorIndex) % config % apply_reg_limits = .true. ! if true application of profiles limits
544 tvs_opts(sensorIndex) % config % verbose = .false. ! verbose output
545 tvs_opts(sensorIndex) % config % do_checkinput = .true. ! to check if input profiles are within absolute and regression limits
546 tvs_opts(sensorIndex) % config % fix_hgpl = .false. ! for backward compatibility with RTTOV-12 should be changed later
547 !< General RT options
548 tvs_opts(sensorIndex) % rt_all % switchrad = .true. ! to use brightness temperature (true) or radiance (false) units in AD routine
549 tvs_opts(sensorIndex) % rt_all % use_q2m = .false. ! if true use of surface humidity (false for compatibility with the way rttov 8.7 was compiled)
550 tvs_opts(sensorIndex) % rt_all % addrefrac = .true. ! to account for atmospheric refraction
551 tvs_opts(sensorIndex) % rt_all % dtau_test = .true. ! for backward compatibility with RTTOV-12 may be changed later
552 tvs_opts(sensorIndex) % rt_all % use_t2m_opdep = .false. ! for backward compatibility with RTTOV-12 may be changed later
553 !< VIS/IR RT options
554 tvs_opts(sensorIndex) % rt_ir % addsolar = .false. ! to model solar component in the near IR (2000 cm-1 et plus)
555 tvs_opts(sensorIndex) % rt_ir % addaerosl = .false. ! to account for scattering due to aerosols
556 tvs_opts(sensorIndex) % rt_ir % addclouds = .false. ! to account for scattering due to clouds
557 tvs_opts(sensorIndex) % rt_ir % ir_sea_emis_model = 2 ! ISEM (ir_sea_emis_model 1) useful for GEORAD
558 ! 2 selects IREMIS which is more sophisticated
559 tvs_opts(sensorIndex) % rt_ir % pc % ipcreg = -1 ! index of the required PC predictors... to see later
560 tvs_opts(sensorIndex) % rt_ir % pc % addpc = .false. ! to carry out principal component calculations
561 tvs_opts(sensorIndex) % rt_ir % pc % addradrec = .false. ! to reconstruct radiances from principal components
562 !< MW RT options
563 tvs_opts(sensorIndex) % rt_mw % clw_data = tvs_isInstrumUsingCLW(tvs_instruments(sensorIndex)) ! disponibilite du profil d'eau liquide
564 tvs_opts(sensorIndex) % rt_mw % fastem_version = 6 ! use fastem version 6 microwave sea surface emissivity model (1-6)
565 tvs_opts(sensorIndex) % rt_mw % clw_scheme = 1 ! default and recommended is 2 just for backward compatibility
566 !< Interpolation options
567 tvs_opts(sensorIndex) % interpolation % addinterp = .true. ! use of internal profile interpolator (rt calculation on model levels)
568 tvs_opts(sensorIndex) % interpolation % lgradp = .true. ! allow tl/ad of user pressure levels
569 tvs_opts(sensorIndex) % interpolation % interp_mode = interp_rochon_loglinear_wfn ! see table 9 page 37 of RTTOV 12.1 users guide
570 tvs_opts(sensorIndex) % interpolation % reg_limit_extrap = tvs_regLimitExtrap
571
572 tvs_opts(sensorIndex) % rt_all % co2_data = .false.
573 tvs_opts(sensorIndex) % rt_all % n2o_data = .false.
574 tvs_opts(sensorIndex) % rt_all % co_data = .false.
575 tvs_opts(sensorIndex) % rt_all % ch4_data = .false.
576
577 if (runObsOperatorWithHydrometeors) then
578 tvs_opts_scatt(sensorIndex) % interp_mode = tvs_opts(sensorIndex) % interpolation % interp_mode ! Set interpolation method
579 tvs_opts_scatt(sensorIndex) % reg_limit_extrap = tvs_regLimitExtrap
580 tvs_opts_scatt(sensorIndex) % fastem_version = tvs_opts(sensorIndex) % rt_mw % fastem_version
581 tvs_opts_scatt(sensorIndex) % supply_foam_fraction = .false.
582 tvs_opts_scatt(sensorIndex) % use_t2m_opdep = tvs_opts(sensorIndex) % rt_all % use_t2m_opdep
583 tvs_opts_scatt(sensorIndex) % use_q2m = tvs_opts(sensorIndex) % rt_all % use_q2m
584 tvs_opts_scatt(sensorIndex) % lgradp = .true.
585 tvs_opts_scatt(sensorIndex) % lusercfrac = .false. !< Switch to enable user-specified effective cloud fraction ??
586 tvs_opts_scatt(sensorIndex) % config % do_checkinput = tvs_opts(sensorIndex) % config % do_checkinput
587 tvs_opts_scatt(sensorIndex) % config % apply_reg_limits = tvs_opts(sensorIndex) % config % apply_reg_limits
588 tvs_opts_scatt(sensorIndex) % config % verbose = .true.
589 tvs_opts_scatt(sensorIndex) % config % fix_hgpl= tvs_opts(sensorIndex) % config % fix_hgpl
590 ! other option may be considered:
591 !real(jprb) :: cc_threshold = 1.E-3_jprb !< Minimum effective cloud fraction threshold to consider scattering
592 !real(jprb) :: ice_polarisation = 1.40_jprb !< Polarised scattering factor for ice hydrometeors (<0 = no polarisation)
593 !logical(jplm) :: ozone_data = .false. !< Switch to enable input of O3 profile
594 ! because standard RTTOV coefficients in the MW have no ozone sensitivity
595 !logical(jplm) :: rad_down_lin_tau = .true. !< Linear-in-tau or layer-mean for downwelling radiances
596 !logical(jplm) :: hydro_cfrac_tlad = .true. !< Switch for hydrometeor TL/AD sensitivity to effective cfrac
597 !logical(jplm) :: zero_hydro_tlad = .false. !< Switch for hydrometeor TL/AD sensitivity in layers with zero
598 ! hydrometeor concentration
599 end if
600
601
602 errorStatus = errorStatus_success
603 call utl_tmg_start(16,'----RttovSetup')
604 write(*,*) ' sensorIndex,tvs_nchan(sensorIndex)', sensorIndex,tvs_nchan(sensorIndex)
605 if ( tvs_mpiTask0ReadCoeffs ) then
606 call tvs_rttov_read_coefs(errorStatus, tvs_coefs(sensorIndex), tvs_opts(sensorIndex), &
607 tvs_ichan(1:tvs_nchan(sensorIndex),sensorIndex), tvs_listSensors(:,sensorIndex))
608 else
609 call rttov_read_coefs ( &
610 errorStatus, &! out
611 tvs_coefs(sensorIndex), &
612 tvs_opts(sensorIndex), &
613 instrument= tvs_listSensors(:,sensorIndex), &! in
614 channels= tvs_ichan(1:tvs_nchan(sensorIndex),sensorIndex) ) ! in option
615 end if
616 if (errorStatus /= errorStatus_success) then
617 write(*,*) 'rttov_read_coefs: fatal error reading coefficients',errorStatus,sensorIndex,tvs_listSensors(1:3,sensorIndex)
618 call utl_abort('tvs_setupAlloc')
619 end if
620
621 if (runObsOperatorWithHydrometeors) then
622 hydrotableFilename = 'hydrotable_' // trim(platform_name(tvs_platforms(sensorIndex))) // '_' // &
623 trim(inst_name(tvs_instruments(sensorIndex))) // '.dat'
624 call rttov_read_scattcoeffs(errorStatus, tvs_opts_scatt(sensorIndex), tvs_coefs(sensorIndex), &
625 tvs_coef_scatt(sensorIndex), file_coef=hydrotableFilename)
626 if (errorStatus /= errorStatus_success) then
627 write(*,*) 'rttov_read_scattcoeffs: fatal error reading RTTOV-SCATT coefficients', hydrotableFilename
628 call utl_abort('tvs_setupAlloc')
629 end if
630 end if
631 call utl_tmg_stop(16)
632
633 tvs_opts(sensorIndex) % rt_all % ozone_data = ( tvs_coefs(sensorIndex) % coef % nozone > 0 ) ! profil d'ozone disponible
634
635 ! Ensure the options and coefficients are consistent
636 call rttov_user_options_checkinput(errorStatus, tvs_opts(sensorIndex), tvs_coefs(sensorIndex))
637 if (errorStatus /= errorStatus_success) then
638 write(*,*) 'error in rttov options', errorStatus
639 call utl_abort('tvs_setupAlloc')
640 end if
641
642 end do
643
644
645 ! 4. Memory allocations for radiative tranfer model variables
646
647 ! Radiance by profile
648
649 allocate( tvs_radiance(tvs_nobtov) ,stat= allocStatus(1))
650
651 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_setupAlloc radiances 1')
652
653 do tovsIndex = 1, tvs_nobtov
654 sensorIndex = tvs_lsensor(tovsIndex)
655 if (sensorIndex > -1) then
656 ! allocate BT equivalent to total direct, tl and ad radiance output
657 allocate( tvs_radiance(tovsIndex) % bt ( tvs_nchan(sensorIndex) ) ,stat= allocStatus(1))
658 tvs_radiance(tovsIndex) % bt ( : ) = 0.d0
659 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_setupAlloc radiances 2')
660 nullify (tvs_radiance(tovsIndex) % clear )
661 end if
662 end do
663
664 end if
665
666 write(*,*) 'Leaving tvs_setupAlloc'
667
668 end subroutine tvs_setupAlloc
669
670 !--------------------------------------------------------------------------
671 ! tvs_getProfile
672 !--------------------------------------------------------------------------
673 subroutine tvs_getProfile(profiles, profileType, cld_profiles_opt)
674 !
675 !:Purpose: sets profiles as a pointer of type rttov_profile
676 ! based on profileType equal to nl or tlad.
677 !
678 implicit none
679
680 ! Arguments:
681 type(rttov_profile), pointer, intent(inout) :: profiles(:)
682 character(len=*), intent(in) :: profileType
683 type(rttov_profile_cloud), pointer, optional, intent(inout) :: cld_profiles_opt(:)
684
685 select case( trim( profileType) )
686 case('nl')
687 profiles => tvs_profiles_nl
688 if (present(cld_profiles_opt)) cld_profiles_opt => tvs_cld_profiles_nl
689 case('tlad')
690 profiles => tvs_profiles_tlad
691 if (present(cld_profiles_opt)) cld_profiles_opt => tvs_cld_profiles_tlad
692 case default
693 call utl_abort('tvs_getProfile: invalid profileType ' // profileType )
694 end select
695
696 end subroutine tvs_getProfile
697
698 !--------------------------------------------------------------------------
699 ! tvs_allocTransmission
700 !--------------------------------------------------------------------------
701 subroutine tvs_allocTransmission(nlevels)
702 !
703 !:Purpose: Allocate the global rttov transmission structure used
704 ! when this is needed for some purpose (e.g. used in
705 ! LETKF to determine peak pressure level of each radiance
706 ! channel for vertical localization).
707 !
708 implicit none
709
710 ! Arguments:
711 integer, intent(in) :: nlevels
712
713 ! Locals:
714 integer :: allocStatus(2), jo, isens, nc
715
716 allocStatus(:) = 0
717 allocate( tvs_transmission(tvs_nobtov), stat=allocStatus(1))
718 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_allocTransmission')
719
720 do jo = 1, tvs_nobtov
721 isens = tvs_lsensor(jo)
722 nc = tvs_nchan(isens)
723 ! allocate transmittance from surface and from pressure levels
724 allocate( tvs_transmission(jo) % tau_total(nc), stat= allocStatus(1))
725 allocate( tvs_transmission(jo) % tau_levels(nlevels,nc), stat= allocStatus(2))
726 call utl_checkAllocationStatus(allocStatus, ' tvs_allocTransmission')
727 end do
728
729 end subroutine tvs_allocTransmission
730
731
732
733 !--------------------------------------------------------------------------
734 ! tvs_setup
735 !--------------------------------------------------------------------------
736 subroutine tvs_setup
737 !
738 !:Purpose: to read namelist NAMTOV, initialize the observation error covariance and setup RTTOV-12.
739 !
740 implicit none
741
742 ! Locals:
743 integer sensorIndex, ierr, nulnam
744 integer, external :: fclos, fnom
745 integer :: instrumentIndex, numMWInstrumToUseCLW, numMWInstrumToUseHydrometeors
746
747 ! Namelist variables: (local)
748 character(len=8) :: crtmodl ! For now, must equal 'RTTOV'
749 integer :: nsensors ! MUST NOT BE INCLUDED IN NAMELIST!
750 character(len=15) :: csatid(tvs_maxNumberOfSensors) ! List of satellite names
751 character(len=15) :: cinstrumentid(tvs_maxNumberOfSensors) ! List of incstrument names
752 logical :: ldbgtov ! Choose to print simulated and observed Tb to listing
753 logical :: useO3Climatology ! Choose to use ozone climatology (otherwise model field)
754 logical :: regLimitExtrap ! Choose to use RTTOV reg_limit_extrap option
755 logical :: doAzimuthCorrection(tvs_maxNumberOfSensors) ! Choose to apply correction to azimuth angle
756 logical :: isAzimuthValid(tvs_maxNumberOfSensors) ! Indicate if azimuth angle is valid
757 logical :: userDefinedDoAzimuthCorrection ! Indicate if user defined azimuth correction is to be used
758 logical :: userDefinedIsAzimuthValid ! Indicate if user defined azimuth angle is valid
759 logical :: mpiTask0ReadCoeffs ! Choose to read coeffs only on task 0 and broadcast
760 logical :: mwInstrumUsingCLW_tl ! Choose to use CLW increment in TL/AD if exists as state variable
761 logical :: mwInstrumUsingHydrometeors_tl ! Choose to all hydomet variables in TL/AD if exist as state variables
762 real(8) :: cloudScaleFactor ! Scale factor applied to model produced clouds to account for bias
763 character(len=15) :: instrumentNamesUsingCLW(tvs_maxNumberOfSensors) ! List of inst names using CLW
764 character(len=15) :: instrumentNamesUsingHydrometeors(tvs_maxNumberOfSensors) ! List of inst name using full set of hydromet variables
765 logical :: mwAllskyAssim ! High-level key to activate all-sky treatment of MW radiances
766
767 namelist /NAMTOV/ nsensors, csatid, cinstrumentid
768 namelist /NAMTOV/ ldbgtov,useO3Climatology
769 namelist /NAMTOV/ useUofWIREmiss, crtmodl
770 namelist /NAMTOV/ useMWEmissivityAtlas, mWAtlasId
771 namelist /NAMTOV/ mwInstrumUsingCLW_tl, instrumentNamesUsingCLW
772 namelist /NAMTOV/ mwInstrumUsingHydrometeors_tl, instrumentNamesUsingHydrometeors
773 namelist /NAMTOV/ regLimitExtrap, doAzimuthCorrection, userDefinedDoAzimuthCorrection
774 namelist /NAMTOV/ isAzimuthValid, userDefinedIsAzimuthValid, cloudScaleFactor
775 namelist /NAMTOV/ mwAllskyAssim, mpiTask0ReadCoeffs
776
777 ! return if the NAMTOV does not exist
778 if ( .not. utl_isNamelistPresent('NAMTOV','./flnml') ) then
779 write(*,*)
780 write(*,*) 'tvs_setup: Namelist block NAMTOV is missing in the namelist.'
781 write(*,*) ' Skipping tvs_setup.'
782 return
783 end if
784
785 ! 1.1 Default values for namelist variables
786
787 nsensors = MPC_missingValue_INT
788 csatid(:) = '***UNDEFINED***'
789 cinstrumentid(:) = '***UNDEFINED***'
790 doAzimuthCorrection(:) = .false.
791 isAzimuthValid(:) = .false.
792 ldbgtov = .false.
793 useO3Climatology = .true.
794 userDefinedDoAzimuthCorrection = .false.
795 userDefinedIsAzimuthValid = .false.
796 crtmodl = 'RTTOV'
797 useUofWIREmiss = .false.
798 useMWEmissivityAtlas = .false.
799 mWAtlasId = 1 !Default to TELSEM-2
800 mwInstrumUsingCLW_tl = .false.
801 mwInstrumUsingHydrometeors_tl = .false.
802 instrumentNamesUsingCLW(:) = '***UNDEFINED***'
803 instrumentNamesUsingHydrometeors(:) = '***UNDEFINED***'
804 regLimitExtrap = .true.
805 cloudScaleFactor = 0.5D0
806 mwAllskyAssim = .false.
807 mpiTask0ReadCoeffs = .true.
808
809 ! 1.2 Read the NAMELIST NAMTOV to modify them
810
811 nulnam = 0
812 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
813 read(nulnam, nml=namtov, iostat=ierr)
814 if (ierr /= 0) call utl_abort('tvs_setup: Error reading namelist NAMTOV')
815
816 if (mmpi_myid == 0) write(*,nml=namtov)
817 ierr = fclos(nulnam)
818
819 ! 1.3 Transfer namelist variables to module variables
820 if (nsensors /= MPC_missingValue_INT) then
821 call utl_abort('tvs_setup: check namelist section NAMTOV; nsensors should be removed as it is' // &
822 ' now computed by Midas from cinstrumentid and csatid arrays')
823 end if
824
825 tvs_nsensors = 0
826 sensor_loop:do sensorIndex = 1, tvs_maxNumberOfSensors
827 if (cinstrumentid(sensorIndex) /= "***UNDEFINED***" .and. &
828 csatid(sensorIndex) /= "***UNDEFINED***" ) then
829 tvs_nsensors = tvs_nsensors + 1
830 else
831 exit sensor_loop
832 end if
833 end do sensor_loop
834
835 tvs_debug = ldbgtov
836 radiativeTransferCode = crtmodl
837 tvs_useO3Climatology = useO3Climatology
838 tvs_instrumentName(:) = cinstrumentid(:)
839 tvs_satelliteName(:) = csatid(:)
840 tvs_mwInstrumUsingCLW_tl = mwInstrumUsingCLW_tl
841 tvs_regLimitExtrap = regLimitExtrap
842 tvs_userDefinedDoAzimuthCorrection = userDefinedDoAzimuthCorrection
843 tvs_userDefinedIsAzimuthValid = userDefinedIsAzimuthValid
844 tvs_doAzimuthCorrection(:) = doAzimuthCorrection(:)
845 tvs_isAzimuthValid(:) = isAzimuthValid(:)
846 tvs_cloudScaleFactor = cloudScaleFactor
847 tvs_mwAllskyAssim = mwAllskyAssim
848 tvs_mpiTask0ReadCoeffs = mpiTask0ReadCoeffs
849
850 ! 1.4 Validate namelist values
851
852 if ( tvs_nsensors == 0 ) then
853 if(mmpi_myid==0) then
854 write(*,*) ' ====================================================='
855 write(*,*) ' tvs_setup: Number of sensors is zero, skipping setup'
856 write(*,*) ' ====================================================='
857 end if
858 return
859 end if
860
861 if ( radiativeTransferCode /= 'RTTOV' ) then
862 write(*,'(A)') ' Invalid radiation model name'
863 call utl_abort('tvs_setup')
864 end if
865
866 if ( tvs_nsensors > tvs_maxNumberOfSensors ) then
867 write(*,'(A)') ' Number of sensors (tvs_nsensors) is greater than maximum allowed (tvs_maxNumberOfSensors)'
868 call utl_abort('tvs_setup')
869 end if
870
871 ! 1.5 Print the content of this NAMELIST
872
873 if(mmpi_myid == 0) then
874 write(*,'(A)')
875 write(*,'(3X,A)') '- Parameters used for TOVS processing (read in NAMTOV)'
876 write(*,'(3X,A)') ' ----------------------------------------------------'
877 write(*,'(6X,A,2X,L1)') 'TOVS debug : ', tvs_debug
878 write(*,'(6X,A,2X,L1)') 'Use of UW IR land emissivity atlases : ', useUofWIREmiss
879 write(*,'(6X,A,2X,L1)') 'Use of MW land emissivity atlases : ', useMWEmissivityAtlas
880 if (useMWEmissivityAtlas) then
881 write(*,'(6X,A,2X,I1)') 'MW atlas Id : ', mWAtlasId
882 end if
883 write(*,'(6X,A,2X,L1)') 'Use of reg_limit_extrap : ', regLimitExtrap
884 write(*,'(6X,A,2X,A)') 'Radiative transfer model : ', radiativeTransferCode
885 write(*,'(6X,A,2X,I3)') 'Number of sensors : ', tvs_nsensors
886 write(*,"(6X,'Satellite ids : ',10A10)") (tvs_satelliteName(sensorIndex), sensorIndex=1,tvs_nsensors)
887 write(*,"(6X,'Instrument ids : ',10A10)") (tvs_instrumentName(sensorIndex), sensorIndex=1,tvs_nsensors)
888 write(*,'(A)')
889 write(*,'(A)')
890 write(*,'(A)')
891 write(*,'(3X,A)') '- Reading and initialization in preparation to the TOVS processing'
892 write(*,'(5X,A)') '----------------------------------------------------------------'
893 end if
894
895 ! 1.6 Set up platform, satellite, instrument and channel mapping
896
897 call sensors()
898
899 ! Get the name and number of instruments to use CLW
900 instrumentIdsUsingCLW(:) = -1
901 do instrumentIndex = 1, tvs_nsensors
902 instrumentIdsUsingCLW(instrumentIndex) = tvs_getInstrumentId(instrumentNamesUsingCLW(instrumentIndex))
903 if ( instrumentNamesUsingCLW(instrumentIndex) /= '***UNDEFINED***' ) then
904 if ( instrumentIdsUsingCLW(instrumentIndex) == -1 ) then
905 write(*,*) instrumentIndex, instrumentNamesUsingCLW(instrumentIndex)
906 call utl_abort('tvs_setup: Unknown instrument name to use CLW')
907 end if
908 else
909 numMWInstrumToUseCLW = instrumentIndex - 1
910 exit
911 end if
912 end do
913
914 ! Get the name and number of instruments to use hydrometeors
915 instrumentIdsUsingHydrometeors(:) = -1
916 do instrumentIndex = 1, tvs_nsensors
917 instrumentIdsUsingHydrometeors(instrumentIndex) = &
918 tvs_getInstrumentId(instrumentNamesUsingHydrometeors(instrumentIndex))
919 if ( instrumentNamesUsingHydrometeors(instrumentIndex) /= '***UNDEFINED***' ) then
920 if ( instrumentIdsUsingHydrometeors(instrumentIndex) == -1 ) then
921 write(*,*) instrumentIndex, instrumentNamesUsingHydrometeors(instrumentIndex)
922 call utl_abort('tvs_setup: Unknown instrument name to use hydrometeors')
923 end if
924 else
925 numMWInstrumToUseHydrometeors = instrumentIndex - 1
926 exit
927 end if
928 end do
929
930 ! check instrument is either using CLW or hydrometeors for non-ATMS instruments
931 do instrumentIndex = 1, numMWInstrumToUseHydrometeors
932 if (numMWInstrumToUseCLW == 0) exit
933
934 if (any(instrumentIdsUsingCLW(1:numMWInstrumToUseCLW) == &
935 instrumentIdsUsingHydrometeors(instrumentIndex))) then
936 write(*,*) instrumentIndex, instrumentNamesUsingHydrometeors(instrumentIndex)
937 call utl_abort('tvs_setup: all-sky TtHu for this intrument is not assimilated yet')
938 end if
939 end do
940
941 tvs_numMWInstrumUsingCLW = numMWInstrumToUseCLW
942 tvs_numMWInstrumUsingHydrometeors = numMWInstrumToUsehydrometeors
943
944 if ( mmpi_myid == 0 ) then
945 write(*,*) 'tvs_setup: Instrument IDs to use CLW: ', instrumentIdsUsingCLW(1:numMWInstrumToUseCLW)
946 write(*,*) 'tvs_setup: Number of instruments to use CLW: ', numMWInstrumToUseCLW
947
948 write(*,*) 'tvs_setup: Instrument IDs to use hydrometeors: ', &
949 instrumentIdsUsingHydrometeors(1:numMWInstrumToUseHydrometeors)
950 write(*,*) 'tvs_setup: Number of instruments to use hydrometeors: ', &
951 numMWInstrumToUseHydrometeors
952 end if
953
954 end subroutine tvs_setup
955
956 !--------------------------------------------------------------------------
957 ! tvs_cleanup
958 !--------------------------------------------------------------------------
959 subroutine tvs_cleanup
960 !
961 !:Purpose: release memory used by RTTOV-12.
962 !
963 implicit none
964
965 ! Locals:
966 integer :: allocStatus(8)
967 integer :: iSensor,iObs,nChans,nl
968
969 Write(*,*) 'tvs_cleanup: Starting'
970
971 allocStatus(:) = 0
972
973 if ( radiativeTransferCode == 'RTTOV' ) then
974
975 !___ radiance by profile
976
977 do iObs = 1, tvs_nobtov
978 iSensor = tvs_lsensor(iObs)
979 nchans = tvs_nchan(isensor)
980 ! deallocate BT equivalent to total direct, tl and ad radiance output
981 deallocate( tvs_radiance(iObs) % bt ,stat= allocStatus(1))
982 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_cleanup radiances 1',.false.)
983 end do
984
985 deallocate( tvs_radiance ,stat= allocStatus(1))
986 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_cleanup radiances 2')
987
988 do iObs = 1, tvs_nobtov
989 iSensor = tvs_lsensor(iObs)
990 nl = tvs_coefs(iSensor) % coef % nlevels
991 ! deallocate model profiles atmospheric arrays with RTTOV levels dimension
992 call rttov_alloc_prof(allocStatus(1),1,tvs_profiles_nl(iObs),nl, & ! 1 = nprofiles un profil a la fois
993 tvs_opts(iSensor),asw=0,coefs=tvs_coefs(iSensor),init=.false. ) ! asw =0 deallocation
994 call rttov_alloc_prof(allocStatus(2),1,tvs_profiles_tlad(iObs),nl, & ! 1 = nprofiles un profil a la fois
995 tvs_opts(iSensor),asw=0,coefs=tvs_coefs(iSensor),init=.false. ) ! asw =0 deallocation
996 call utl_checkAllocationStatus(allocStatus(1:2), 'profiles deallocation in tvs_cleanup',.false.)
997 end do
998
999 deallocate(tvs_profiles_nl, stat=allocStatus(1) )
1000 deallocate(tvs_profiles_tlad, stat=allocStatus(2) )
1001 call utl_checkAllocationStatus(allocStatus(1:2), ' tvs_setupAlloc tvs_profiles_nl/tlad')
1002
1003 do iSensor = tvs_nsensors,1,-1
1004 call rttov_dealloc_coefs(allocStatus(1), tvs_coefs(iSensor) )
1005 Write(*,*) 'Deallocating coefficient structure for instrument', iSensor
1006 call utl_checkAllocationStatus(allocStatus(1:1), ' rttov_dealloc_coefs in tvs_cleanup', .false.)
1007 end do
1008
1009 deallocate (tvs_coefs ,stat= allocStatus(1))
1010 deallocate (tvs_listSensors ,stat= allocStatus(2))
1011 deallocate (tvs_opts ,stat= allocStatus(3))
1012
1013 call utl_checkAllocationStatus(allocStatus(1:3), ' tvs_cleanup', .false.)
1014
1015 end if
1016
1017 deallocate (tvs_nchan, stat= allocStatus(1))
1018 deallocate (tvs_ichan, stat= allocStatus(2))
1019 deallocate (tvs_lsensor, stat= allocStatus(3))
1020 deallocate (tvs_headerIndex, stat= allocStatus(4))
1021 deallocate (tvs_tovsIndex, stat= allocStatus(5))
1022 deallocate (tvs_isReallyPresent,stat= allocStatus(6))
1023 deallocate (tvs_nchanMpiGlobal, stat= allocStatus(7))
1024 deallocate (tvs_ichanMpiGlobal, stat= allocStatus(8))
1025
1026 call utl_checkAllocationStatus(allocStatus, ' tvs_cleanup', .false.)
1027
1028 Write(*,*) 'tvs_cleanup: Finished'
1029
1030 end subroutine tvs_cleanup
1031
1032 !--------------------------------------------------------------------------
1033 ! tvs_deallocateProfilesNlTlAd
1034 !--------------------------------------------------------------------------
1035 subroutine tvs_deallocateProfilesNlTlAd
1036 !
1037 !:Purpose: release memory used by RTTOV-12.
1038 !
1039 implicit none
1040
1041 ! Locals:
1042 integer :: allocStatus(8)
1043
1044 Write(*,*) 'tvs_deallocateProfilesNlTlAd: Starting'
1045
1046 allocStatus(:) = 0
1047
1048 if ( radiativeTransferCode == 'RTTOV' ) then
1049 if ( allocated(tvs_profiles_nl) ) deallocate(tvs_profiles_nl, stat=allocStatus(1))
1050 if ( allocated(tvs_profiles_tlad) ) deallocate(tvs_profiles_tlad, stat=allocStatus(2))
1051 call utl_checkAllocationStatus(allocStatus(1:2), ' tvs_profiles_nl tvs_profiles_tlad', .false.)
1052 end if
1053
1054 Write(*,*) 'tvs_deallocateProfilesNlTlAd: Finished'
1055
1056 end subroutine tvs_deallocateProfilesNlTlAd
1057
1058 !--------------------------------------------------------------------------
1059 ! sensors
1060 !--------------------------------------------------------------------------
1061 subroutine sensors
1062 !
1063 !:Purpose: Initialisation of the RTTOV-10 platform, satellite
1064 ! and instrument ID's. Also set burp to RTTOV channel
1065 ! mapping offset.
1066 ! To verify and transfom the sensor information contained in the
1067 ! NAMTOV namelist into the variables required by RTTTOV-7:
1068 ! platform, satellite and instrument ID's.
1069 !
1070 implicit none
1071
1072 ! Locals:
1073 integer sensorIndex, instrumentIndex, platformIndex
1074 integer ipos1, ipos2
1075 integer numerosat, ierr, kindex, nulnam
1076 character(len=15) :: tempocsatid
1077 logical, save :: first=.true.
1078 integer, save :: ioffset1b(0:ninst-1)
1079 character(len=15) :: tempo_inst
1080 integer, external :: fnom, fclos
1081
1082 ! Namelist variables:
1083 character(len=8) :: listinstrum(0:ninst-1) ! List of instrument names
1084 integer :: listoffset(0:ninst-1) ! Corresponding list of channel offset values
1085 namelist /NAMCHANOFFSET/ listoffset, listinstrum
1086
1087 ! 1.0 Go through sensors and set RTTOV-10 variables
1088
1089 do sensorIndex=1, tvs_nsensors
1090 tvs_platforms (sensorIndex) = -1
1091 tvs_satellites (sensorIndex) = -1
1092 tvs_instruments(sensorIndex) = -1
1093 tvs_channelOffset(sensorIndex) = -1
1094 end do
1095
1096 if (first) then
1097 if ( utl_isNamelistPresent('NAMCHANOFFSET', './flnml') ) then
1098 call utl_abort('sensors: NAMCHANOFFSET namelist section should be now in flnml_static !')
1099 end if
1100 ! read the namelist
1101 nulnam = 0
1102 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1103 if (ierr /= 0) then
1104 write(*,*) 'Error while opening flnml_static namelist file !'
1105 call utl_abort('sensors')
1106 end if
1107 listoffset(:) = 0
1108 listinstrum(:) = 'XXXXXXXX'
1109 read(nulnam,NAMCHANOFFSET, iostat=ierr)
1110 if (ierr /= 0) then
1111 write(*,*) 'Error while reading NAMCHANOFFSET namelist section in flnml_static file !'
1112 call utl_abort('sensors')
1113 end if
1114 do instrumentIndex=0, ninst - 1
1115 if ( listinstrum(instrumentIndex) /= 'XXXXXXXX' ) then
1116 ioffset1b( tvs_getInstrumentId( listinstrum(instrumentIndex) ) ) = listoffset(instrumentIndex)
1117 end if
1118 end do
1119 ierr = fclos(nulnam)
1120 first = .false.
1121 end if
1122
1123 ! 1.1 Set platforms and satellites
1124
1125 ! N.B.: Special cases for satellites TERRA and AQUA.
1126 ! For consistency with the RTTOV-10 nomenclature, rename:
1127 ! TERRA to eos1
1128 ! AQUA to eos2
1129 ! NPP to jpss0
1130 ! JPSS to jpss0
1131 ! HMWARI to himawari
1132 ! FY-3C to FY3-3
1133 do sensorIndex = 1, tvs_nsensors
1134 if ( tvs_satelliteName(sensorIndex) == 'TERRA' ) then
1135 tempocsatid = 'eos1'
1136 else if ( tvs_satelliteName(sensorIndex) == 'AQUA' ) then
1137 tempocsatid = 'eos2'
1138 else if ( tvs_satelliteName(sensorIndex) == 'NPP' ) then
1139 tempocsatid = 'jpss0'
1140 else if ( tvs_satelliteName(sensorIndex) == 'JPSS' ) then
1141 tempocsatid = 'jpss0'
1142 else if ( tvs_satelliteName(sensorIndex)(1:6) == 'HMWARI' ) then
1143 tempocsatid = 'himawari' // trim(tvs_satelliteName(sensorIndex) (7:15))
1144 else if ( tvs_satelliteName(sensorIndex) == 'FY-3C' ) then
1145 TEMPOCSATID = 'FY3-3'
1146 else
1147 call up2low(tvs_satelliteName(sensorIndex),tempocsatid)
1148 end if
1149 do platformIndex = 1, nplatforms
1150 ipos1 = len_trim(platform_name(platformIndex))
1151 ipos2 = index(tempocsatid,platform_name(platformIndex)(1:ipos1))
1152 if ( ipos2 == 1 ) then
1153 tvs_platforms(sensorIndex) = platformIndex
1154 kindex = platformIndex
1155 end if
1156 end do
1157 if ( tvs_platforms(sensorIndex) < 0 ) then
1158 write(*,'(A)') ' Satellite ' // trim(tempocsatid) // ' not supported.'
1159 call utl_abort('SENSORS')
1160 else
1161 ipos1 = len_trim(platform_name(kindex))
1162 ipos2 = len_trim(tempocsatid)
1163 read(tempocsatid(ipos1+1:ipos2),*,IOSTAT=ierr) numerosat
1164 numerosat = abs ( numerosat )
1165 if ( ierr /= 0) then
1166 write(*,'(A,1x,i6,1x,i3,1x,i3,1x,A15)') 'Problem while reading satellite number', &
1167 ierr, ipos1, ipos2, tempocsatid
1168 call utl_abort('SENSORS')
1169 else
1170 tvs_satellites(sensorIndex) = numerosat
1171 end if
1172 end if
1173 end do
1174
1175 ! 1.2 Set instruments,
1176 ! also set channel offset, which is in fact a channel mapping between
1177 ! the channel number in BURP files and the channel number used in
1178 ! RTTOV-10.
1179
1180 do sensorIndex = 1, tvs_nsensors
1181 if ( tvs_instrumentName(sensorIndex)(1:10) == 'GOESIMAGER') then !cas particulier
1182 tvs_instruments(sensorIndex) = inst_id_goesim
1183 else if ( tvs_satelliteName(sensorIndex)(1:5) == 'MTSAT') then !autre cas particulier
1184 tvs_instruments(sensorIndex) = inst_id_gmsim
1185 else
1186 call up2low(tvs_instrumentName(sensorIndex),tempo_inst)
1187 do instrumentIndex = 0, ninst -1
1188 if ( trim(tempo_inst) == trim(inst_name(instrumentIndex))) then
1189 tvs_instruments(sensorIndex) = instrumentIndex
1190 end if
1191 end do
1192 end if
1193 if ( tvs_instruments(sensorIndex) < 0 ) then
1194 write(*,'(A)') ' INSTRUMENT '// trim( tvs_instrumentName(sensorIndex)) // ' not supported.'
1195 call utl_abort('SENSORS')
1196 end if
1197 tvs_channelOffset(sensorIndex) = ioffset1b(tvs_instruments(sensorIndex))
1198 end do
1199
1200 ! 1.3 Print the RTTOV-12 related variables
1201
1202 if (mmpi_myid == 0) then
1203 write(*,*)
1204 write(*,'(3X,A)') '- SENSORS. Variables prepared for RTTOV-12:'
1205 write(*,'(3X,A)') ' ----------------------------------------'
1206 write(*,*)
1207 write(*,'(6X,A,I3)') 'Number of sensors : ', tvs_nsensors
1208 write(*,"('Platform numbers : ',6X,10I3)") (tvs_platforms(sensorIndex), sensorIndex=1,tvs_nsensors)
1209 write(*,"('Satellite numbers : ',6X,10I3)") (tvs_satellites(sensorIndex), sensorIndex=1,tvs_nsensors)
1210 write(*,"('Instrument numbers : ',6X,10I3)") (tvs_instruments(sensorIndex), sensorIndex=1,tvs_nsensors)
1211 write(*,"('Channel mapping offsets : ',6X,10I3)") (tvs_channelOffset(sensorIndex), sensorIndex=1,tvs_nsensors)
1212 end if
1213
1214 end subroutine sensors
1215
1216
1217 !--------------------------------------------------------------------------
1218 ! tvs_getAllIdBurpTovs
1219 !--------------------------------------------------------------------------
1220 subroutine tvs_getAllIdBurpTovs(idatypListSize, idatypList)
1221 !
1222 !:Purpose: Function to return a list of all idatyp (a.k.a. codtyp) values
1223 ! for all possible radiance observations (according to the namelist)
1224 !
1225 implicit none
1226
1227 ! Argument:
1228 integer, intent(out) :: idatypListSize
1229 integer, intent(out) :: idatypList(:)
1230
1231 ! Locals:
1232 logical, save :: first=.true.
1233 integer, save :: ninst_tovs
1234 integer :: nulnam, ierr, instrumentIndex
1235 integer, external :: fnom, fclos
1236 integer, save :: list_inst(ninst)
1237
1238 ! Namelist variables:
1239 character(len=22) :: inst_names(ninst) ! List of instrument names for all radiance types
1240 namelist /namtovsinst/ inst_names
1241
1242 if (tvs_nsensors == 0) then
1243 ! no tovs data will be read, therefore false
1244 idatypList(:) = MPC_missingValue_int
1245 idatypListSize = 0
1246 return
1247 end if
1248
1249 if (first) then
1250 if ( utl_isNamelistPresent('NAMTOVSINST', './flnml') ) then
1251 call utl_abort('tvs_getAllIdBurpTovs: NAMTOVSINST namelist section should be now in flnml_static !')
1252 end if
1253 nulnam = 0
1254 ninst_tovs = 0
1255 list_inst(:) = -1
1256 inst_names(:) = 'XXXXXX'
1257 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1258 read(nulnam, nml=namtovsinst, iostat=ierr)
1259 if (ierr /= 0) call utl_abort('tvs_getAllIdBurpTovs: Error reading NAMTOVSINST namelist section in flnml_static file')
1260 if (mmpi_myid == 0) write(*,nml=namtovsinst)
1261 ierr = fclos(nulnam)
1262 do instrumentIndex=1, ninst
1263 if (inst_names(instrumentIndex) == 'XXXXXX' ) then
1264 ninst_tovs = instrumentIndex - 1
1265 exit
1266 else
1267 list_inst(instrumentIndex) = codtyp_get_codtyp( inst_names(instrumentIndex) )
1268 if (list_inst(instrumentIndex) < 0) then
1269 write(*,*) inst_names(instrumentIndex)
1270 call utl_abort('tvs_isIdBurpTovs: unknown instrument in namtovsinst namelist')
1271 end if
1272 end if
1273 end do
1274 if ( ninst_tovs == 0 ) call utl_abort('tvs_getAllIdBurpTovs: Empty namtovsinst namelist')
1275 first = .false.
1276 end if
1277
1278 idatypList(:) = MPC_missingValue_int
1279 idatypList(1:ninst_tovs) = list_inst(1:ninst_tovs)
1280 idatypListSize = ninst_tovs
1281
1282 end subroutine tvs_getAllIdBurpTovs
1283
1284 !--------------------------------------------------------------------------
1285 ! tvs_isIdBurpTovs
1286 !--------------------------------------------------------------------------
1287 logical function tvs_isIdBurpTovs(idatyp)
1288 !
1289 !:Purpose: Function to check if the given idatyp (a.k.a. codtyp)
1290 ! corresponds to a radiance
1291 !
1292 implicit none
1293
1294 ! Arguments:
1295 integer, intent(in) :: idatyp
1296
1297 ! Locals:
1298 logical, save :: first=.true.
1299 integer, save :: ninst_tovs
1300 integer :: nulnam, ierr, instrumentIndex
1301 integer, external :: fnom, fclos
1302 integer, save :: list_inst(ninst)
1303
1304 ! Namelist variables:
1305 character(len=22) :: inst_names(ninst) ! List of instrument names for all radiance types
1306 namelist /namtovsinst/ inst_names
1307
1308 if (tvs_nsensors == 0) then
1309 ! no tovs data will be read, therefore false
1310 tvs_isIdBurpTovs = .false.
1311 return
1312 end if
1313
1314 if (first) then
1315 if ( utl_isNamelistPresent('NAMTOVSINST', './flnml') ) then
1316 call utl_abort('tvs_isIdBurpTovs: NAMTOVSINST namelist section should be now in flnml_static !')
1317 end if
1318 nulnam = 0
1319 ninst_tovs = 0
1320 list_inst(:) = -1
1321 inst_names(:) = 'XXXXXX'
1322 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1323 read(nulnam, nml=namtovsinst, iostat=ierr)
1324 if (ierr /= 0) call utl_abort('tvs_isIdBurpTovs: Error reading NAMTOVSINST namelist section in flnml_static file')
1325 if (mmpi_myid == 0) write(*,nml=namtovsinst)
1326 ierr = fclos(nulnam)
1327 do instrumentIndex=1, ninst
1328 if (inst_names(instrumentIndex) == 'XXXXXX' ) then
1329 ninst_tovs= instrumentIndex - 1
1330 exit
1331 else
1332 list_inst(instrumentIndex) = codtyp_get_codtyp( inst_names(instrumentIndex) )
1333 if (list_inst(instrumentIndex) < 0) then
1334 write(*,*) inst_names(instrumentIndex)
1335 call utl_abort('tvs_isIdBurpTovs: unknown instrument in namtovsinst namelist')
1336 end if
1337 end if
1338 end do
1339 if ( ninst_tovs == 0 ) call utl_abort('tvs_isIdBurpTovs: Empty namtovsinst namelist')
1340 first = .false.
1341 end if
1342
1343 tvs_isIdBurpTovs = .false.
1344
1345 do instrumentIndex = 1, ninst_tovs
1346 if (idatyp == list_inst(instrumentIndex) ) then
1347 tvs_isIdBurpTovs = .true.
1348 exit
1349 end if
1350 end do
1351
1352 end function tvs_isIdBurpTovs
1353
1354 !--------------------------------------------------------------------------
1355 ! tvs_isIdBurpHyperSpectral
1356 !--------------------------------------------------------------------------
1357 logical function tvs_isIdBurpHyperSpectral(idatyp)
1358 !
1359 !:Purpose: Function to check if the given idatyp (a.k.a. codtyp)
1360 ! corresponds to a hyper-spectral infrared radiance
1361 !
1362 implicit none
1363
1364 ! Argument:
1365 integer, intent(in) :: idatyp
1366
1367 ! Locals:
1368 logical, save :: first=.true.
1369 integer, save :: ninst_hyper
1370 integer :: nulnam, ierr, instrumentIndex
1371 integer, external :: fnom, fclos
1372 integer, save :: list_inst(ninst)
1373
1374 ! Namelist variables:
1375 character(len=22) :: name_inst(ninst) ! List of instrument names for hyperspectral IR
1376 namelist /namhyper/ name_inst
1377
1378 if (tvs_nsensors == 0) then
1379 ! no tovs data will be read, therefore false
1380 tvs_isIdBurpHyperSpectral = .false.
1381 return
1382 end if
1383
1384 if (first) then
1385 if ( utl_isNamelistPresent('NAMHYPER', './flnml') ) then
1386 call utl_abort('tvs_isIdBurpHyperSpectral: NAMHYPER namelist section should be now in flnml_static !')
1387 end if
1388 nulnam = 0
1389 ninst_hyper = 0
1390 list_inst(:) = -1
1391 name_inst(:) = 'XXXXXX'
1392 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1393 read(nulnam, nml=namhyper, iostat=ierr)
1394 if (ierr /= 0) call utl_abort('tvs_isIdBurpHyperSpectral: Error reading NAMHYPER namelist section in flnml_static file')
1395 if (mmpi_myid == 0) write(*,nml=namhyper)
1396 ierr = fclos(nulnam)
1397 do instrumentIndex=1, ninst
1398 if (name_inst(instrumentIndex) == 'XXXXXX' ) then
1399 ninst_hyper = instrumentIndex - 1
1400 exit
1401 else
1402 list_inst(instrumentIndex) = codtyp_get_codtyp( name_inst(instrumentIndex) )
1403 if (list_inst(instrumentIndex) < 0) then
1404 write(*,*) name_inst(instrumentIndex)
1405 call utl_abort('tvs_isIdBurpHyperSpectral: unknown instrument in namhyper namelist')
1406 end if
1407 end if
1408 end do
1409 if ( ninst_hyper == 0 ) call utl_abort('tvs_isIdBurpHyperSpectral: Empty namhyper namelist')
1410 first = .false.
1411 end if
1412
1413 tvs_isIdBurpHyperSpectral = .false.
1414
1415 do instrumentIndex = 1, ninst_hyper
1416 if (idatyp == list_inst(instrumentIndex) ) then
1417 tvs_isIdBurpHyperSpectral = .true.
1418 exit
1419 end if
1420 end do
1421
1422 end function tvs_isIdBurpHyperSpectral
1423
1424 !--------------------------------------------------------------------------
1425 ! tvs_isIdBurpInst
1426 !--------------------------------------------------------------------------
1427 logical function tvs_isIdBurpInst(idburp,cinst)
1428 !
1429 !:Purpose: function to check if the provided idburp (a.k.a. codtyp) corresponds to instrument cinst
1430 !
1431 implicit none
1432
1433 ! Arguments:
1434 integer, intent(in) :: idburp ! Input codtyp
1435 character(len=*), intent(in) :: cinst ! Input instrument name
1436
1437 if (tvs_nsensors == 0) then
1438 ! no tovs data will be read, therefore false
1439 tvs_isIdBurpInst = .false.
1440 return
1441 end if
1442
1443 tvs_isIdBurpInst = ( idburp == codtyp_get_codtyp(cinst) )
1444
1445 end function tvs_isIdBurpInst
1446
1447 !--------------------------------------------------------------------------
1448 ! tvs_getPlatformId
1449 !--------------------------------------------------------------------------
1450 integer function tvs_getPlatformId(name)
1451 !
1452 !:Purpose: return RTTOV platform id (>0) from platform name.
1453 ! -1 if not found
1454 !
1455 implicit none
1456
1457 ! Arguments:
1458 character(len=*), intent(in) :: name ! Platform name
1459
1460 ! Locals:
1461 integer :: platformIndex, length, ipos
1462 character(len=64) :: tempo_name
1463
1464 tvs_getPlatformId = -1
1465 length = len_trim(name)
1466 call up2low(name(1:length),tempo_name(1:length))
1467
1468 if ( index(tempo_name(1:length),'npp') /= 0 ) then
1469 tvs_getPlatformId = platform_id_jpss
1470 else if ( index(tempo_name(1:length),'hmwari') /= 0 ) then
1471 tvs_getPlatformId = platform_id_himawari
1472 else
1473 do platformIndex = 1, nplatforms
1474 ipos = index(tempo_name(1:length),trim(platform_name(platformIndex)))
1475 if (ipos == 1) then
1476 tvs_getPlatformId = platformIndex
1477 exit
1478 end if
1479 end do
1480 end if
1481
1482 end function tvs_getPlatformId
1483
1484 !--------------------------------------------------------------------------
1485 ! tvs_getInstrumentId
1486 !--------------------------------------------------------------------------
1487 integer function tvs_getInstrumentId(name)
1488 !
1489 !:Purpose: return RTTOV instrument id from intrument name. 0 is a valid answer.
1490 ! -1 if not found
1491 !
1492 implicit none
1493
1494 ! Arguments:
1495 character(len=*), intent(in) :: name ! Instrument name
1496
1497 ! Locals:
1498 integer :: instrumentIndex, length
1499 character(len=64) :: tempo_name
1500
1501 tvs_getInstrumentId = -1
1502 length = len_trim(name)
1503 call up2low(name(1:length),tempo_name(1:length))
1504 if ( trim(tempo_name(1:length)) == 'goesim' ) then
1505 tvs_getInstrumentId = inst_id_goesim
1506 else if ( trim(tempo_name(1:length)) == 'gmsim' ) then
1507 tvs_getInstrumentId = inst_id_gmsim
1508 else if ( trim(tempo_name(1:length)) == 'mtsatim' ) then
1509 tvs_getInstrumentId = inst_id_mtsatim
1510 else
1511 do instrumentIndex = 0, ninst - 1
1512 if (trim(inst_name(instrumentIndex)) == tempo_name(1:length) ) then
1513 tvs_getInstrumentId = instrumentIndex
1514 exit
1515 end if
1516 end do
1517 end if
1518 end function tvs_getInstrumentId
1519
1520 !--------------------------------------------------------------------------
1521 ! tvs_isInstrumHyperSpectral
1522 !--------------------------------------------------------------------------
1523 logical function tvs_isInstrumHyperSpectral(instrum)
1524 !
1525 !:Purpose: given an RTTOV instrument code return if it is an hyperspectral one
1526 ! information from namelist NAMHYPER
1527 !
1528 implicit none
1529
1530 ! Arguments:
1531 integer, intent(in) :: instrum ! input Rttov instrument code
1532
1533 ! Locals:
1534 integer, parameter :: maxsize = 100
1535 integer :: nulnam, ierr, instrumentIndex
1536 integer, save :: list_inst(maxsize), ninst_hir
1537 logical, save :: first = .true.
1538 integer, external :: fclos, fnom
1539
1540 ! Namelist variables:
1541 character (len=8) :: name_inst(maxsize) ! List of instrument names for hyperspectral IR
1542 namelist /NAMHYPER/ name_inst
1543
1544 if (first) then
1545 if ( utl_isNamelistPresent('NAMHYPER', './flnml') ) then
1546 call utl_abort('tvs_isInstrumHyperSpectral: NAMHYPER namelist section should be now in flnml_static !')
1547 end if
1548 nulnam = 0
1549 ninst_hir = 0
1550 name_inst(:) = 'XXXXXXX'
1551 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1552 read(nulnam,nml=namhyper, iostat=ierr)
1553 if (ierr /= 0) call utl_abort('tvs_isInstrumHyperSpectral: Error reading namelist section NAMHYPER in flnm_static file')
1554 if (mmpi_myid == 0) write(*,nml=namhyper)
1555 ierr = fclos(nulnam)
1556 list_inst(:) = -1
1557 do instrumentIndex=1, maxsize
1558 list_inst(instrumentIndex) = tvs_getInstrumentId( name_inst(instrumentIndex) )
1559 if (name_inst(instrumentIndex) /= 'XXXXXXX') then
1560 if (list_inst(instrumentIndex) == -1) then
1561 write(*,*) instrumentIndex,name_inst(instrumentIndex)
1562 call utl_abort('tvs_isInstrumHyperSpectral: Unknown instrument name')
1563 end if
1564 else
1565 ninst_hir = instrumentIndex - 1
1566 exit
1567 end if
1568 end do
1569 first = .false.
1570 if (ninst_hir == 0) then
1571 write(*,*) 'tvs_isInstrumHyperSpectral: Warning : empty namhyper namelist !'
1572 end if
1573 end if
1574 tvs_isInstrumHyperSpectral = .false.
1575 do instrumentIndex =1, ninst_hir
1576 if ( instrum == list_inst(instrumentIndex)) then
1577 tvs_isInstrumHyperSpectral = .true.
1578 exit
1579 end if
1580 end do
1581
1582 end function tvs_isInstrumHyperSpectral
1583
1584 !--------------------------------------------------------------------------
1585 ! tvs_isNameHyperSpectral
1586 !--------------------------------------------------------------------------
1587 logical function tvs_isNameHyperSpectral(cinstrum)
1588 !
1589 !:Purpose: given an instrument name
1590 ! returns if it is an hyperspectral one
1591 ! (information from namelist NAMHYPER)
1592 !
1593 implicit none
1594
1595 ! Arguments:
1596 character(len=*), intent(in) :: cinstrum
1597
1598 ! Locals:
1599 integer, parameter :: maxsize = 20
1600 integer :: nulnam, ierr, i
1601 integer, save :: ninst_hir
1602 logical, save :: lfirst = .true.
1603 integer, external :: fclos, fnom
1604 character (len=8) :: name2
1605
1606 ! Namelist variables:
1607 character (len=8),save :: name_inst(maxsize) ! List of instrument names for hyperspectral IR
1608 namelist /NAMHYPER/ name_inst
1609
1610 if (lfirst) then
1611 nulnam = 0
1612 ninst_hir = 0
1613 name_inst(:) = 'XXXXXXX'
1614 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1615 read(nulnam,nml=namhyper, iostat=ierr)
1616 if (ierr /= 0) call utl_abort('tvs_isNameHyperSpectral: Error reading NAMHYPER namelist section in flnml_static file')
1617 if (mmpi_myid == 0) write(*,nml=namhyper)
1618 ierr = fclos(nulnam)
1619 do i=1, maxsize
1620 if (name_inst(i) == 'XXXXXXX') then
1621 ninst_hir = i -1
1622 exit
1623 end if
1624 end do
1625 lfirst = .false.
1626 if (ninst_hir == 0) then
1627 write(*,*) 'tvs_isNameHyperSpectral: Warning : empty namhyper namelist !'
1628 end if
1629 end if
1630
1631 tvs_isNameHyperSpectral = .false.
1632
1633 call up2low(cinstrum, name2)
1634
1635 do i=1, ninst_hir
1636 if ( trim(name2) == name_inst(i)) then
1637 tvs_isNameHyperSpectral = .true.
1638 exit
1639 end if
1640 end do
1641
1642 end function tvs_isNameHyperSpectral
1643
1644 !--------------------------------------------------------------------------
1645 ! tvs_isInstrumGeostationary
1646 !--------------------------------------------------------------------------
1647 logical function tvs_isInstrumGeostationary(instrum)
1648 !
1649 !:Purpose: given an RTTOV instrument code return if it is a Geostationnary Imager
1650 ! information from namelist NAMGEO
1651 !
1652 implicit none
1653
1654 ! Arguments:
1655 integer, intent(in) :: instrum ! input Rttov instrument code
1656
1657 ! Locals:
1658 integer, parameter :: maxsize = 100
1659 integer :: nulnam, ierr, instrumentIndex
1660 integer, save :: list_inst(maxsize), ninst_geo
1661 logical, save :: first = .true.
1662 integer, external :: fnom, fclos
1663
1664 ! Namelist variables:
1665 character(len=8) :: name_inst(maxsize) ! List of instrument names for geostationary
1666 namelist /NAMGEO/ name_inst
1667
1668 if (first) then
1669 if ( utl_isNamelistPresent('NAMGEO', './flnml') ) then
1670 call utl_abort('tvs_isInstrumGeostationary: NAMGEO namelist section should be now in flnml_static !')
1671 end if
1672 nulnam = 0
1673 ninst_geo = 0
1674 name_inst(:) = 'XXXXXX'
1675 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1676 read(nulnam,nml=namgeo, iostat=ierr)
1677 if (ierr /= 0) call utl_abort('tvs_isInstrumGeostationary: Error reading namelist section NAMGEO in flnml_static file')
1678 if (mmpi_myid == 0) write(*,nml=namgeo)
1679 ierr = fclos(nulnam)
1680 list_inst(:) = -1
1681 do instrumentIndex=1, maxsize
1682 list_inst(instrumentIndex) = tvs_getInstrumentId( name_inst(instrumentIndex) )
1683 if (name_inst(instrumentIndex) /= 'XXXXXX') then
1684 if (list_inst(instrumentIndex) == -1) then
1685 write(*,*) instrumentIndex,name_inst(instrumentIndex)
1686 call utl_abort('tvs_isInstrumGeostationary: Unknown instrument name')
1687 end if
1688 else
1689 ninst_geo = instrumentIndex - 1
1690 exit
1691 end if
1692 end do
1693 first = .false.
1694 if (ninst_geo == 0) then
1695 write(*,*) 'tvs_isInstrumGeostationary: Warning : empty namgeo namelist !'
1696 end if
1697 end if
1698 tvs_isInstrumGeostationary = .false.
1699 do instrumentIndex=1, ninst_geo
1700 if ( instrum == list_inst(instrumentIndex)) then
1701 tvs_isInstrumGeostationary = .true.
1702 exit
1703 end if
1704 end do
1705
1706 end function tvs_isInstrumGeostationary
1707
1708 !--------------------------------------------------------------------------
1709 ! tvs_isInstrumUsingCLW
1710 !--------------------------------------------------------------------------
1711 function tvs_isInstrumUsingCLW(instrumId) result(idExist)
1712 !
1713 !:Purpose: given an RTTOV instrument code return if it is in the list to use CLW
1714 !
1715 implicit none
1716
1717 ! Arguments:
1718 integer, intent(in) :: instrumId ! input Rttov instrument code
1719 ! Result:
1720 logical :: idExist
1721
1722 ! Locals:
1723 integer :: instrumentIndex
1724
1725 idExist = .false.
1726 do instrumentIndex = 1, tvs_numMWInstrumUsingCLW
1727 if ( instrumId == instrumentIdsUsingCLW(instrumentIndex) ) then
1728 idExist = .true.
1729 exit
1730 end if
1731 end do
1732
1733 end function tvs_isInstrumUsingCLW
1734
1735 !--------------------------------------------------------------------------
1736 ! tvs_isInstrumUsingHydrometeors
1737 !--------------------------------------------------------------------------
1738 function tvs_isInstrumUsingHydrometeors(instrumId) result(idExist)
1739 !
1740 !:Purpose: given an RTTOV instrument code return if it is in the list to use Hydrometeors
1741 !
1742 implicit none
1743
1744 ! Arguments:
1745 integer, intent(in) :: instrumId ! input Rttov instrument code
1746 ! Result:
1747 logical :: idExist
1748
1749 ! Locals:
1750 integer :: instrumentIndex
1751
1752 idExist = .false.
1753 do instrumentIndex = 1, tvs_numMWInstrumUsingHydrometeors
1754 if ( instrumId == instrumentIdsUsingHydrometeors(instrumentIndex) ) then
1755 idExist = .true.
1756 exit
1757 end if
1758 end do
1759
1760 end function tvs_isInstrumUsingHydrometeors
1761
1762 !--------------------------------------------------------------------------
1763 ! tvs_isInstrumAllskyTtAssim
1764 !--------------------------------------------------------------------------
1765 function tvs_isInstrumAllskyTtAssim(instrumId) result(allskyTtAssim)
1766 !
1767 !:Purpose: determine if all-sky temperature-channel assimilation is asked for the instrument.
1768 !
1769 implicit none
1770
1771 ! Arguments:
1772 integer, intent(in) :: instrumId ! input Rttov instrument code
1773 ! Result:
1774 logical :: allskyTtAssim
1775
1776 allskyTtAssim = (tvs_mwAllskyAssim .and. tvs_isInstrumUsingCLW(instrumId))
1777
1778 end function tvs_isInstrumAllskyTtAssim
1779
1780 !--------------------------------------------------------------------------
1781 ! tvs_isInstrumAllskyHuAssim
1782 !--------------------------------------------------------------------------
1783 function tvs_isInstrumAllskyHuAssim(instrumId) result(allskyHuAssim)
1784 !
1785 !:Purpose: determine if all-sky humidity-channel assimilation is asked for the instrument.
1786 !
1787 implicit none
1788
1789 ! Arguments:
1790 integer, intent(in) :: instrumId ! input Rttov instrument code
1791 ! Result:
1792 logical :: allskyHuAssim
1793
1794 allskyHuAssim = (tvs_mwAllskyAssim .and. tvs_isInstrumUsingHydrometeors(instrumId))
1795
1796 end function tvs_isInstrumAllskyHuAssim
1797
1798 !--------------------------------------------------------------------------
1799 ! tvs_mapInstrum
1800 !--------------------------------------------------------------------------
1801 subroutine tvs_mapInstrum(instrumburp,instrum)
1802 !
1803 !:Purpose: Map burp satellite instrument (element #2019) to RTTOV-7 instrument.
1804 ! A negative value is returned, if no match in found.
1805 !
1806 !:Table of RTTOV-7 instrument identifier:
1807 !
1808 ! ================== ===================== ==================
1809 ! Instrument Instrument identifier Sensor type
1810 ! ================== ===================== ==================
1811 ! HIRS 0 ir
1812 ! MSU 1 mw
1813 ! SSU 2 ir
1814 ! AMSUA 3 mw
1815 ! AMSUB 4 mw
1816 ! AVHRR 5 ir
1817 ! SSMI 6 mw
1818 ! VTPR1 7 ir
1819 ! VTPR2 8 ir
1820 ! TMI 9 mw
1821 ! SSMIS 10 mw
1822 ! AIRS 11 ir
1823 ! MODIS 13 ir
1824 ! ATSR 14 ir
1825 ! MHS 15 mw
1826 ! ATMS 19 mw
1827 ! MVIRI 20 ir
1828 ! SEVIRI 21 ir
1829 ! GOESIMAGER 22 ir
1830 ! GOESSOUNDER 23 ir
1831 ! GMS/MTSAT IMAGER 24 ir
1832 ! FY2-VISSR 25 ir
1833 ! FY1-MVISR 26 ir
1834 ! AHI 56 ir
1835 ! ================== ===================== ==================
1836 !
1837 implicit none
1838
1839 ! Arguments:
1840 integer, intent(in) :: instrumburp ! burp satellite instrument (element #2019)
1841 integer, intent(out) :: instrum ! RTTOV-7 instrument ID numbers (e.g. 3 for AMSUA)
1842
1843 ! Locals:
1844 integer instrumentIndex, numinstburp
1845 integer, parameter :: mxinstrumburp = 100
1846 logical, save :: first = .true.
1847 integer :: nulnam, ier
1848 integer, external :: fnom, fclos
1849
1850 ! Namelist variables:
1851 integer, save :: listburp(mxinstrumburp) ! List of instrument ID values from obs file
1852 character(len=8), save :: listinstrum(mxinstrumburp) ! List of instrument names
1853 namelist /NAMINST/ listburp, listinstrum
1854
1855 ! Table of BURP satellite sensor identifier element #002019
1856
1857 ! 1.0 Find instrument
1858
1859 if (first) then
1860 if ( utl_isNamelistPresent('NAMINST', './flnml') ) then
1861 call utl_abort('tvs_mapInstrum: NAMINST namelist section should be now in flnml_static !')
1862 end if
1863
1864 ! set the default values
1865 listburp(:) = -1
1866 listinstrum(:) = 'XXXXXXXX'
1867
1868 ! read the namelist
1869 nulnam = 0
1870 ier = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1871 if (ier /= 0) then
1872 write(*,*) 'Error while opening flnml_static namelist file !'
1873 call utl_abort('tvs_mapInstrum')
1874 end if
1875 read(nulnam,NAMINST,iostat=ier)
1876 if (ier /= 0) then
1877 write(*,*) 'Error while reading NAMINST namelist section in flnml_static file !'
1878 call utl_abort('tvs_mapInstrum')
1879 end if
1880 ier = fclos(nulnam)
1881
1882 ! figure out how many valid elements in the lists
1883 do instrumentIndex=1, mxinstrumburp
1884 if(listburp(instrumentIndex) == -1) then
1885 numinstburp = instrumentIndex - 1
1886 exit
1887 end if
1888 end do
1889 if (numinstburp > mxinstrumburp) then
1890 call utl_abort('tvs_mapInstrum: exceeded maximum number of platforms')
1891 end if
1892 write(*,*) 'tvs_mapInstrum: number of satellites found in namelist = ',numinstburp
1893 write(*,*) 'tvs_mapInstrum: listburp = ',listburp(1:numinstburp)
1894 write(*,*) 'tvs_mapInstrum: listinstrum = ',listinstrum(1:numinstburp)
1895 first=.false.
1896 end if
1897
1898 instrum = -1
1899 do instrumentIndex=1, mxinstrumburp
1900 if ( instrumburp == listburp(instrumentIndex) ) then
1901 instrum = tvs_getInstrumentId( listinstrum(instrumentIndex) )
1902 exit
1903 end if
1904 end do
1905
1906 end subroutine tvs_mapInstrum
1907
1908 !--------------------------------------------------------------------------
1909 ! tvs_isNameGeostationary
1910 !--------------------------------------------------------------------------
1911 logical function tvs_isNameGeostationary(cinstrum)
1912 !
1913 !:Purpose: given an instrument name following BUFR convention
1914 ! returns if it is a Geostationnary Imager
1915 ! (information from namelist NAMGEOBUFR)
1916 !
1917 implicit none
1918
1919 ! Arguments:
1920 character(len=*), intent(in) :: cinstrum
1921
1922 ! Locals:
1923 integer, parameter :: maxsize = 100
1924 integer :: nulnam, ierr, i
1925 integer, save :: ninst_geo
1926 logical, save :: lfirst = .true.
1927 integer, external :: fnom, fclos
1928
1929 ! Namelist variables:
1930 character (len=8),save :: name_inst(maxsize) ! List of instrument names for geostationary
1931 namelist /NAMGEOBUFR/ name_inst
1932
1933 if (lfirst) then
1934 if ( utl_isNamelistPresent('NAMGEOBUFR', './flnml') ) then
1935 call utl_abort('tvs_isNameGeostationary: NAMGEOBUFR namelist section should be now in flnml_static !')
1936 end if
1937 nulnam = 0
1938 ninst_geo = 0
1939 name_inst(:) = 'XXXXXXXX'
1940 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
1941 read(nulnam,nml=namgeobufr, iostat=ierr)
1942 if (ierr /= 0) call utl_abort('tvs_isNameGeostationary: Error reading namelist section NAMGEOBUFR in flnml_static_file')
1943 if (mmpi_myid == 0) write(*,nml=namgeobufr)
1944 ierr = fclos(nulnam)
1945 do i=1, maxsize
1946 if (name_inst(i) == 'XXXXXXXX') then
1947 ninst_geo = i - 1
1948 exit
1949 end if
1950 end do
1951 lfirst = .false.
1952 if (ninst_geo == 0) then
1953 write(*,*) 'tvs_isNameGeostationary: Warning : empty namgeobufr namelist !'
1954 end if
1955 end if
1956
1957 tvs_isNameGeostationary = .false.
1958 do i=1, ninst_geo
1959 if ( trim(cinstrum) == trim(name_inst(i)) ) then
1960 tvs_isNameGeostationary= .true.
1961 exit
1962 end if
1963 end do
1964
1965 end function tvs_isNameGeostationary
1966
1967 !--------------------------------------------------------------------------
1968 ! tvs_mapSat
1969 !--------------------------------------------------------------------------
1970 subroutine tvs_mapSat(isatBURP,iplatform,isat)
1971 !
1972 !:Purpose: Map burp satellite identifier (element #1007)
1973 ! to RTTOV-7 platform and satellite.
1974 ! Negative values are returned, if no match in found.
1975 !
1976 !:Table of RTTOV-7 platform identifier:
1977 !
1978 ! ======== ===========================
1979 ! Platform RTTOV-7 platform identifier
1980 ! ======== ===========================
1981 ! NOAA 1
1982 ! DMSP 2
1983 ! METEOSAT 3
1984 ! GOES 4
1985 ! GMS 5
1986 ! FY2 6
1987 ! TRMM 7
1988 ! ERS 8
1989 ! EOS 9
1990 ! METOP 10
1991 ! ENVISAT 11
1992 ! MSG 12
1993 ! FY1 13
1994 ! ADEOS 14
1995 ! MTSAT 15
1996 ! CORIOLIS 16
1997 ! NPP 17
1998 ! ======== ===========================
1999 !
2000 !:Example:
2001 ! NOAA15, which has a burp satellite identifier value of 206,
2002 ! is mapped into the following:
2003 ! RTTOV-7 platform = 1,
2004 ! RTTOV-7 satellite = 15.
2005 !
2006 !:Arguments:
2007 ! :isatBURP: BURP satellite identifier
2008 ! :iplatform: RTTOV-7 platform ID numbers (e.g. 1 for NOAA)
2009 ! :isat: RTTOV-7 satellite ID numbers (e.g. 15)
2010 !
2011 implicit none
2012
2013 ! Arguments:
2014 integer, intent(in) :: isatburp ! BURP satellite identifier
2015 integer, intent(out) :: iplatform ! RTTOV-7 platform ID numbers (e.g. 1 for NOAA)
2016 integer, intent(out) :: isat ! RTTOV-7 satellite ID numbers (e.g. 15)
2017
2018 ! Locals:
2019 integer :: satelliteIndex, ierr, nulnam
2020 logical, save :: first=.true.
2021 integer, external :: fnom, fclos
2022 integer, parameter:: mxsatburp = 100
2023 integer, save :: numsatburp
2024
2025 ! Namelist variables:
2026 integer, save :: listburp(mxsatburp) ! Table of BURP satellite identifier element #001007
2027 character(len=8), save :: listplat(mxsatburp) ! Table of RTTOV platform identifier
2028 integer, save :: listsat (mxsatburp) ! Table of RTTOV satellite identifier
2029
2030 namelist /NAMSAT/ listburp, listplat, listsat
2031
2032 ! Fill tables from namelist at the first call
2033 if (first) then
2034 if ( utl_isNamelistPresent('NAMSAT', './flnml') ) then
2035 call utl_abort('tvs_mapSat: NAMSAT namelist section should be now in flnml_static !')
2036 end if
2037 ! set the default values
2038 listburp(:) = -1
2039 listsat(:) = -1
2040 listplat(:) = 'XXXXXXXX'
2041 ! read the namelist
2042 nulnam = 0
2043 ierr = fnom(nulnam,'./flnml_static','FTN+SEQ+R/O',0)
2044 if (ierr /= 0) then
2045 write(*,*) 'Error while opening namelist flnml_static file !'
2046 call utl_abort('tvs_mapSat')
2047 end if
2048 read(nulnam, NAMSAT, iostat = ierr)
2049 if (ierr /= 0) then
2050 write(*,*) 'Error while reading NAMSAT namelist section in flnml_static file !'
2051 call utl_abort('tvs_mapSat')
2052 end if
2053 ierr = fclos(nulnam)
2054
2055 ! Figure out how many valid elements in the lists
2056 do satelliteIndex=1, mxsatburp
2057 if(listburp(satelliteIndex) == -1) then
2058 numsatburp = satelliteIndex - 1
2059 exit
2060 end if
2061 end do
2062 if(numsatburp >= mxsatburp) then
2063 call utl_abort('tvs_mapSat: exceeded maximum number of platforms')
2064 end if
2065 write(*,*) 'tvs_mapSat: number of satellites found in namelist = ',numsatburp
2066 write(*,*) 'tvs_mapSat: listburp = ',listburp(1:numsatburp)
2067 write(*,*) 'tvs_mapSat: listsat = ',listsat(1:numsatburp)
2068 write(*,*) 'tvs_mapSat: listplat = ',listplat(1:numsatburp)
2069 first = .false.
2070 end if
2071
2072 ! Find platform and satellite
2073
2074 iplatform = -1
2075 isat = -1
2076 do satelliteIndex=1, numsatburp
2077 if ( isatburp == listburp(satelliteIndex) ) then
2078 iplatform = tvs_getPlatformId( listplat(satelliteIndex) )
2079 isat = listsat (satelliteIndex)
2080 exit
2081 end if
2082 end do
2083
2084 end subroutine tvs_mapSat
2085
2086 !--------------------------------------------------------------------------
2087 ! tvs_getChanProf
2088 !--------------------------------------------------------------------------
2089 subroutine tvs_getChanprof(sensorTovsIndexes, obsSpaceData, chanprof, lchannel_subset_opt, iptobs_cma_opt)
2090 !
2091 !:Purpose: subroutine to initialize the chanprof structure used by RTTOV
2092 !
2093 implicit none
2094
2095 ! Arguments:
2096 integer, intent(in) :: sensorTovsIndexes(:)
2097 type(struct_obs), intent(in) :: obsSpaceData
2098 type(rttov_chanprof), intent(out) :: chanprof(:)
2099 logical, optional, intent(out) :: lchannel_subset_opt(:,:)
2100 integer, optional, intent(out) :: iptobs_cma_opt(:)
2101
2102 ! Locals:
2103 integer :: count, profileIndex, headerIndex, istart, iend, bodyIndex, channelNumber, iobs
2104 integer :: ChannelIndex
2105
2106 ! Build the list of channels/profiles indices
2107 count = 0
2108 if (present( lchannel_subset_opt )) lchannel_subset_opt(:,:) = .false.
2109
2110 do profileIndex = 1, size(sensorTovsIndexes)
2111 iobs = sensorTovsIndexes(profileIndex)
2112 headerIndex = tvs_headerIndex(iobs)
2113 if (headerIndex > 0) then
2114 istart = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2115 iend= obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + istart - 1
2116 do bodyIndex = istart, iend
2117 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) then
2118 call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
2119 channelNumber, channelIndex )
2120 if (channelIndex >0) then
2121 count = count + 1
2122 chanprof(count)%prof = profileIndex
2123 chanprof(count)%chan = channelIndex
2124 if (present(iptobs_cma_opt)) iptobs_cma_opt(count) = bodyIndex
2125 if (present( lchannel_subset_opt )) lchannel_subset_opt(profileIndex,channelIndex) = .true.
2126 else
2127 write(*,*) 'tvs_getChanProf: strange channel number',channelNumber
2128 end if
2129 end if
2130 end do
2131 end if
2132 end do
2133
2134 end subroutine tvs_getChanprof
2135
2136 !--------------------------------------------------------------------------
2137 ! tvs_countRadiances
2138 !--------------------------------------------------------------------------
2139 integer function tvs_countRadiances(sensorTovsIndexes, obsSpaceData, assim_flag_val_opt)
2140 !
2141 !:Purpose: to count radiances selected for assimilation
2142 !
2143 implicit none
2144
2145 ! Arguments:
2146 integer, intent(in) :: sensorTovsIndexes(:)
2147 type(struct_obs), intent(in) :: obsSpaceData
2148 integer, optional, intent(in) :: assim_flag_val_opt
2149
2150 ! Locals:
2151 integer :: profileIndex, headerIndex, istart, iend, bodyIndex, iobs, assim_flag_val
2152
2153 if (present(assim_flag_val_opt)) then
2154 assim_flag_val = assim_flag_val_opt
2155 else
2156 assim_flag_val = 1
2157 end if
2158
2159 tvs_countRadiances = 0
2160 do profileIndex = 1, size(sensorTovsIndexes)
2161 iobs = sensorTovsIndexes(profileIndex)
2162 headerIndex = tvs_headerIndex(iobs)
2163 if (headerIndex > 0) then
2164 istart = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2165 iend = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + istart - 1
2166 do bodyIndex = istart, iend
2167 if(obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) tvs_countRadiances = tvs_countRadiances + 1
2168 end do
2169 end if
2170 end do
2171
2172 end function tvs_countRadiances
2173
2174 !--------------------------------------------------------------------------
2175 ! tvs_ChangedStypValue(obsspacedata, headerIndex)
2176 !--------------------------------------------------------------------------
2177 integer function tvs_ChangedStypValue(obsSpaceData, headerIndex)
2178 !
2179 !:Purpose: to obtain new STYP value given observed STYP and TTYP value
2180 !
2181 implicit none
2182
2183 ! Arguments:
2184 integer, intent(in) :: headerIndex
2185 type(struct_obs), intent(in) :: obsSpaceData
2186
2187 ! Locals:
2188 integer :: terrainType
2189 integer :: landSea
2190
2191 terrainType = obs_headElem_i(obsSpaceData,OBS_TTYP,headerIndex)
2192 landSea = obs_headElem_i(obsSpaceData,OBS_STYP,headerIndex)
2193
2194 if ( terrainType == 0 ) then
2195 tvs_ChangedStypValue = surftype_seaice
2196 else
2197 tvs_ChangedStypValue = landSea
2198 end if
2199
2200 end function tvs_ChangedStypValue
2201
2202 !--------------------------------------------------------------------------
2203 ! tvs_getHIREmissivities
2204 !--------------------------------------------------------------------------
2205 subroutine tvs_getHIREmissivities(sensorTovsIndexes, obsSpaceData, surfem)
2206 !
2207 !:Purpose: to get emissivity for Hyperspectral Infrared Sounders (AIRS, IASI, CrIS, ...)
2208 !
2209 implicit none
2210
2211 ! Arguments:
2212 integer, intent(in) :: sensorTovsIndexes(:)
2213 type(struct_obs), intent(in) :: obsSpaceData
2214 real(8), intent(out) :: surfem(:)
2215
2216 ! Locals:
2217 integer :: count, profileIndex, iobs, istart, iend, bodyIndex, headerIndex
2218
2219 count = 0
2220 surfem(:) = 0.98d0
2221 do profileIndex = 1, size(sensorTovsIndexes)
2222 iobs = sensorTovsIndexes(profileIndex)
2223 headerIndex = tvs_headerIndex(iobs)
2224 if (headerIndex > 0 ) then
2225 istart = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2226 iend = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + istart - 1
2227 do bodyIndex = istart, iend
2228 if(obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) then
2229 count = count + 1
2230 surfem ( count ) = obs_bodyElem_r(obsSpaceData,OBS_SEM,bodyIndex)
2231 end if
2232 end do
2233 end if
2234 end do
2235
2236 end subroutine tvs_getHIREmissivities
2237
2238 !--------------------------------------------------------------------------
2239 ! tvs_getOtherEmissivities
2240 !--------------------------------------------------------------------------
2241 subroutine tvs_getOtherEmissivities(chanprof, sensorTovsIndexes, sensorType, instrument, surfem, calcemis)
2242 !
2243 !:Purpose: to get emissivity for microwave sounders and infrared geostationary imagers
2244 !
2245 implicit none
2246
2247 ! Arguments:
2248 type(rttov_chanprof), intent(in) :: chanprof(:)
2249 integer, intent(in) :: sensorTovsIndexes(:)
2250 integer, intent(in) :: sensorType
2251 integer, intent(in) :: instrument
2252 real(8), intent(out) :: surfem(:)
2253 logical, intent(out) :: calcemis(:)
2254
2255 ! Locals:
2256 integer :: radiance_index, profileIndex, iobs, surfaceType
2257
2258 do radiance_index = 1, size(chanprof)
2259 profileIndex = chanprof(radiance_index)%prof
2260 iobs = sensorTovsIndexes(profileIndex)
2261 surfaceType = tvs_profiles_nl(iobs) % skin % surftype
2262 if ( sensorType == sensor_id_mw ) then
2263 if ( surfaceType == surftype_land .or. &
2264 surfaceType == surftype_seaice ) then
2265 calcemis(radiance_index) = .false.
2266 surfem (radiance_index) = 0.75d0
2267 else
2268 calcemis(radiance_index) = .true.
2269 surfem (radiance_index) = 0.d0
2270 end if
2271 else if ( tvs_isInstrumHyperSpectral(instrument) ) then
2272 calcemis(radiance_index) = .false.
2273 else if ( tvs_isInstrumGeostationary(instrument) ) then
2274 calcemis(radiance_index) = .true.
2275 surfem (radiance_index) = 0.d0
2276 else
2277 write(*,*) sensorType,instrument
2278 call utl_abort('tvs_getOtherEmissivities. invalid sensor type or unknown IR instrument')
2279 end if
2280 end do
2281
2282 end subroutine tvs_getOtherEmissivities
2283
2284 !--------------------------------------------------------------------------
2285 ! tvs_fillProfiles
2286 !--------------------------------------------------------------------------
2287 subroutine tvs_fillProfiles(columnTrl, obsSpaceData, datestamp, profileType, beSilent)
2288 !
2289 !:Purpose: to fill in tvs_profiles_nl structure before call to non-linear,
2290 ! tangent-linear or adjoint of RTTOV
2291 !
2292 implicit none
2293
2294 ! Arguments:
2295 type(struct_columnData), intent(in) :: columnTrl ! Column structure
2296 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData structure
2297 integer, intent(in) :: datestamp ! CMC date stamp
2298 character(len=*), intent(in) :: profileType
2299 logical, intent(in) :: beSilent ! To control verbosity
2300
2301 ! Locals:
2302 integer :: instrum, iplatform
2303 integer :: nobmax
2304 integer :: sensorIndex, tovsIndex
2305 integer :: profileCount, headerIndex
2306 integer :: profileIndex, levelIndex
2307 integer :: ilowlvl_M,ilowlvl_T,nlv_M,nlv_T
2308 integer :: Vcode
2309 integer :: ierr,day,month,year,ijour,itime
2310 integer :: allocStatus(13)
2311 integer,external :: omp_get_num_threads
2312 integer,external :: newdate
2313 integer, allocatable :: sensorTovsIndexes(:)
2314 integer, allocatable :: sensorHeaderIndexes(:)
2315 type(struct_vco), pointer :: vco
2316 real(8), allocatable :: pressure (:,:)
2317 real(8), allocatable :: latitudes(:)
2318 real(8), allocatable :: ozone(:,:)
2319 character(len=4) :: ozoneVarName
2320 real(8), allocatable :: clw (:,:)
2321 real(8), allocatable :: ciw (:,:)
2322 real(8), allocatable :: rainflux (:,:)
2323 real(8), allocatable :: snowflux (:,:)
2324 real(8), allocatable :: cloudFraction(:,:)
2325 logical, allocatable :: surfTypeIsWater(:)
2326 logical :: runObsOperatorWithClw
2327 logical :: runObsOperatorWithHydrometeors
2328 type(rttov_profile), pointer :: profiles(:)
2329 type(rttov_profile_cloud), pointer :: cld_profiles(:)
2330 real(8), pointer :: column_ptr(:)
2331 real(8) :: zmax, wind
2332
2333 if ( .not. beSilent ) write(*,*) 'tvs_fillProfiles: Starting'
2334
2335 if (tvs_nobtov == 0) return ! exit if there are no tovs data
2336
2337 if ( tvs_numMWInstrumUsingCLW > 0 .and. .not. col_varExist(columnTrl,'LWCR') ) then
2338 call utl_abort('tvs_fillProfiles: if number of instrument to use CLW greater than zero, ' // &
2339 'the LWCR variable must be included as an analysis variable in NAMSTATE. ')
2340 end if
2341 if (tvs_numMWInstrumUsingHydrometeors > 0) then
2342 if (.not. (col_varExist(columnTrl,'LWCR') .and. col_varExist(columnTrl,'IWCR'))) then
2343 call utl_abort('tvs_fillProfiles: if number of instrument to use hydrometeors greater than zero, ' // &
2344 'the LWCR/IWCR variables must be included as an analysis variable in NAMSTATE. ')
2345 end if
2346 if (.not. beSilent .and. .not. (col_varExist(columnTrl,'RF') .and. col_varExist(columnTrl,'SF') .and. &
2347 col_varExist(columnTrl,'CLDR'))) then
2348 write(*,*) 'tvs_fillProfiles: one of RF/SF/CLDR does not exist in NAMSTATE'
2349 end if
2350 end if
2351
2352 if ( (tvs_numMWInstrumUsingCLW == 0 .and. tvs_numMWInstrumUsingHydrometeors == 0 .and. &
2353 tvs_mwAllskyAssim) .or. &
2354 (tvs_numMWInstrumUsingCLW > 0 .and. tvs_numMWInstrumUsingHydrometeors == 0 .and. &
2355 .not. tvs_mwAllskyAssim) .or. &
2356 (tvs_numMWInstrumUsingCLW == 0 .and. tvs_numMWInstrumUsingHydrometeors > 0 .and. &
2357 .not. tvs_mwAllskyAssim) ) then
2358 call utl_abort('tvs_fillProfiles: number of instrument to use CLW/hydrometeors do not match ' // &
2359 'all-sky namelist variable.')
2360 end if
2361
2362 if (.not. tvs_useO3Climatology .and. .not. col_varExist(columnTrl,'TO3') .and. &
2363 .not. col_varExist(columnTrl,'O3L') ) then
2364 call utl_abort('tvs_fillProfiles: if tvs_useO3Climatology is set to .false. the ozone variable ' // &
2365 'must be included as an analysis variable in NAMSTATE. ')
2366 else if (.not.tvs_useO3Climatology) then
2367 if (col_varExist(columnTrl,'TO3') ) then
2368 ozoneVarName = 'TO3'
2369 else
2370 ozoneVarName = 'O3L'
2371 end if
2372 end if
2373
2374 allocStatus(:) = 0
2375
2376 if ( profileType == 'nl' ) then
2377 if ( .not. allocated( tvs_profiles_nl) ) then
2378 allocate(tvs_profiles_nl(tvs_nobtov) , stat=allocStatus(1) )
2379 if (tvs_numMWInstrumUsingHydrometeors > 0) then
2380 allocate(tvs_cld_profiles_nl(tvs_nobtov) , stat=allocStatus(2))
2381 end if
2382 call utl_checkAllocationStatus(allocStatus(1:2), ' tvs_fillProfiles tvs_profiles_nl')
2383 end if
2384 else if ( profileType == 'tlad' ) then
2385 if ( .not. allocated( tvs_profiles_tlad) ) then
2386 allocate(tvs_profiles_tlad(tvs_nobtov) , stat=allocStatus(1) )
2387 if (tvs_numMWInstrumUsingHydrometeors > 0) then
2388 allocate(tvs_cld_profiles_tlad(tvs_nobtov) , stat=allocStatus(2))
2389 end if
2390 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_fillProfiles tvs_profiles_tlad')
2391 else
2392 return
2393 end if
2394 else
2395 write(*,*) 'Invalid profileType ', profileType
2396 call utl_abort('tvs_fillProfiles')
2397 end if
2398
2399 if ( .not. beSilent ) write(*,*) 'tvs_fillProfiles: profileType is ', profileType
2400
2401 call tvs_getProfile(profiles, profileType, cld_profiles)
2402
2403!
2404! 1. Set index for model's lowest level and model top
2405! . ------------------------------------------------
2406
2407 nlv_M = col_getNumLev(columnTrl,'MM')
2408 nlv_T = col_getNumLev(columnTrl,'TH')
2409
2410 if ( col_getPressure(columnTrl,1,1,'TH') < col_getPressure(columnTrl,nlv_T,1,'TH') ) then
2411 ilowlvl_M = nlv_M
2412 ilowlvl_T = nlv_T
2413 else
2414 ilowlvl_M = 1
2415 ilowlvl_T = 1
2416 end if
2417
2418 vco => col_getVco(columnTrl)
2419 Vcode = vco%Vcode
2420
2421 ierr = newdate(datestamp,ijour,itime,-3)
2422 if (ierr < 0) then
2423 write(*,*) 'Invalid datestamp ',datestamp,ijour,itime,ierr
2424 call utl_abort('tvs_fillProfiles')
2425 end if
2426 year= ijour / 10000
2427 month = mod(ijour / 100,100)
2428 day = mod(ijour,100)
2429
2430 ! 1.2 Read ozone climatology
2431
2432 if (tvs_useO3Climatology) call ozo_read_climatology(datestamp)
2433
2434 ! 2. Fill profiles structure
2435
2436 ! loop over all instruments
2437 sensor_loop: do sensorIndex=1, tvs_nsensors
2438
2439 runObsOperatorWithClw = (col_varExist(columnTrl,'LWCR') .and. tvs_numMWInstrumUsingCLW /= 0 .and. &
2440 tvs_isInstrumUsingCLW(tvs_instruments(sensorIndex)))
2441
2442 runObsOperatorWithHydrometeors = (col_varExist(columnTrl,'LWCR') .and. col_varExist(columnTrl,'IWCR') .and. &
2443 tvs_isInstrumUsingHydrometeors(tvs_instruments(sensorIndex)))
2444
2445 if (runObsOperatorWithClw .and. runObsOperatorWithHydrometeors) then
2446 call utl_abort('tvs_fillProfiles: this instrument is mentioned in using CLW and hydrometeors.')
2447 end if
2448
2449 ! first loop over all obs.
2450 profileCount = 0
2451 bobs1: do tovsIndex = 1, tvs_nobtov
2452 if (tvs_lsensor(tovsIndex) == sensorIndex) then
2453 profileCount = profileCount + 1
2454 NOBMAX = tovsIndex
2455 end if
2456 end do bobs1
2457
2458 if (profileCount == 0) cycle sensor_loop
2459
2460 if (tvs_coefs(sensorIndex)%coef%fmv_model_ver >= 9) then
2461 zmax = zenmaxv9
2462 else
2463 zmax = zenmax
2464 end if
2465
2466 allocStatus(:) = 0
2467 allocate (sensorTovsIndexes(profileCount), stat = allocStatus(1) )
2468 allocate (sensorHeaderIndexes(profileCount), stat = allocStatus(2) )
2469 allocate (latitudes(profileCount), stat = allocStatus(3) )
2470 allocate (ozone(nlv_T,profileCount), stat = allocStatus(4) )
2471 allocate (pressure(nlv_T,profileCount), stat = allocStatus(5) )
2472 if (runObsOperatorWithClw .or. runObsOperatorWithHydrometeors) then
2473 allocate (clw (nlv_T,profileCount),stat= allocStatus(6))
2474 clw(:,:) = qlim_getMinValueCloud('LWCR')
2475 end if
2476 if (runObsOperatorWithHydrometeors) then
2477 allocate (ciw (nlv_T,profileCount),stat= allocStatus(7))
2478 allocate (rainFlux (nlv_T,profileCount),stat= allocStatus(8))
2479 allocate (snowFlux (nlv_T,profileCount),stat= allocStatus(9))
2480 allocate (cloudFraction(nlv_T,profileCount),stat= allocStatus(10))
2481 ciw(:,:) = qlim_getMinValueCloud('IWCR')
2482 rainFlux(:,:) = qlim_getMinValueCloud('RF')
2483 snowFlux(:,:) = qlim_getMinValueCloud('SF')
2484 cloudFraction(:,:) = qlim_getMinValueCloud('CLDR')
2485 end if
2486 allocate (surfTypeIsWater(profileCount),stat= allocStatus(11))
2487 surfTypeIsWater(:) = .false.
2488
2489 call utl_checkAllocationStatus(allocStatus, ' tvs_fillProfiles')
2490
2491 profileCount = 0
2492
2493 ! second loop over all obs.
2494 bobs2: do tovsIndex = 1, NOBMAX
2495 if (tvs_lsensor(tovsIndex) /= sensorIndex) cycle bobs2
2496 profileCount = profileCount + 1
2497 sensorTovsIndexes(profileCount) = tovsIndex
2498 headerIndex = tvs_headerIndex(tovsIndex)
2499 sensorHeaderIndexes(profileCount) = headerIndex
2500
2501 call rttov_alloc_prof( &
2502 allocStatus(1), &
2503 1, & ! 1 = nprofiles un profil a la fois
2504 profiles(tovsIndex:tovsIndex), &
2505 nlv_T, &
2506 tvs_opts(sensorIndex), &
2507 asw=1, & ! asw =1 allocation
2508 coefs=tvs_coefs(sensorIndex), &
2509 init=.true. )
2510 if (runObsOperatorWithHydrometeors) then
2511 call rttov_alloc_scatt_prof( &
2512 allocstatus(2), &
2513 1, &
2514 cld_profiles(tovsIndex:tovsIndex), &
2515 nlv_T, &
2516 nhydro=5, & ! depending on what is defined in the Mie tables
2517 nhydro_frac=1, & ! 1 cloud fraction for all variable or nhydro 1 cloud fraction for each variable
2518 asw=1_jpim, & ! 1 => allocate
2519 init=.true., & ! initialize profiles to zero
2520 flux_conversion=[1,2,0,0,0] ) !flux_conversion input units: 0 (default) => kg/kg,
2521 ! 1,2 => kg/m2/s, optional for rain, snow
2522 end if
2523 call utl_checkAllocationStatus(allocStatus(1:2), ' tvs_setupAlloc tvs_fillProfiles')
2524
2525 ! extract land/sea/sea-ice flag (0=land, 1=sea, 2=sea-ice)
2526 profiles(tovsIndex) % skin % surftype = tvs_ChangedStypValue(obsSpaceData,headerIndex)
2527
2528 ! extract satellite zenith and azimuth angle,
2529 ! sun zenith angle, cloud fraction, latitude and longitude
2530 profiles(tovsIndex) % zenangle = obs_headElem_r(obsSpaceData,OBS_SZA,headerIndex)
2531
2532 call validateRttovVariable(profiles(tovsIndex) % zenangle, "satellite zenith angle", &
2533 obsSpaceData, headerIndex, 0.d0, zmax)
2534
2535 profiles(tovsIndex) % azangle = tvs_getCorrectedSatelliteAzimuth(obsSpaceData, headerIndex)
2536 profiles(tovsIndex) % sunazangle = obs_headElem_r(obsSpaceData,OBS_SAZ,headerIndex) ! necessaire pour radiation solaire
2537 iplatform = tvs_coefs(sensorIndex) % coef % id_platform
2538 instrum = tvs_coefs(sensorIndex) % coef % id_inst
2539 profiles(tovsIndex) % sunzenangle = obs_headElem_r(obsSpaceData,OBS_SUN,headerIndex)
2540 if (tvs_opts(sensorIndex)%rt_ir%addsolar) then
2541 call validateRttovVariable(profiles(tovsIndex) % sunzenangle, "sun zenith angle", &
2542 obsSpaceData, headerIndex, 0.d0)
2543 end if
2544 latitudes(profileCount) = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex) * MPC_DEGREES_PER_RADIAN_R8
2545 !1d-5 was chosen as a threshold because it is the high precision for latitudes in BUFR/BURP files
2546 if (latitudes(profileCount) > 90.d0 .and. (latitudes(profileCount)-90.d0) < 1.d-5) latitudes(profileCount) = 90.d0
2547 if (latitudes(profileCount) < -90.d0 .and. (-latitudes(profileCount)-90.d0) < 1.d-5) latitudes(profileCount) = -90.d0
2548 call validateRttovVariable(latitudes(profileCount), 'latitude', obsSpaceData, headerIndex, -90.d0, 90.d0)
2549
2550 profiles(tovsIndex) % longitude = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex) * MPC_DEGREES_PER_RADIAN_R8
2551
2552 surfTypeIsWater(profileCount) = ( tvs_ChangedStypValue(obsSpaceData,headerIndex) == surftype_sea )
2553
2554 do levelIndex = 1, nlv_T
2555 pressure(levelIndex,profileCount) = col_getPressure(columnTrl,levelIndex,headerIndex,'TH') * MPC_MBAR_PER_PA_R8
2556 if ((runObsOperatorWithClw .and. surfTypeIsWater(profileCount)) .or. &
2557 (runObsOperatorWithHydrometeors .and. surfTypeIsWater(profileCount))) then
2558
2559 ! cloud liquid water content
2560 clw(levelIndex,profileCount) = col_getElem(columnTrl,levelIndex,headerIndex,'LWCR')
2561 if (clw(levelIndex,profileCount) < qlim_getMinValueCloud('LWCR') .or. &
2562 clw(levelIndex,profileCount) > qlim_getMaxValueCloud('LWCR')) then
2563 write(*,*) 'tvs_fillProfiles: clw=' , clw(:,profileCount)
2564 call utl_abort('tvs_fillProfiles: columnTrl has clw outside RTTOV bounds')
2565 end if
2566
2567 clw(levelIndex,profileCount) = clw(levelIndex,profileCount) * tvs_cloudScaleFactor
2568 end if
2569 if (runObsOperatorWithHydrometeors .and. surfTypeIsWater(profileCount)) then
2570 ! cloud ice water content
2571 ciw(levelIndex,profileCount) = col_getElem(columnTrl,levelIndex,headerIndex,'IWCR')
2572 if (ciw(levelIndex,profileCount) < qlim_getMinValueCloud('IWCR') .or. &
2573 ciw(levelIndex,profileCount) > qlim_getMaxValueCloud('IWCR')) then
2574 write(*,*) 'tvs_fillProfiles: ciw=' , ciw(:,profileCount)
2575 call utl_abort('tvs_fillProfiles: columnTrl has ciw outside RTTOV bounds')
2576 end if
2577
2578 ! rain flux (zero, if not part of control variables)
2579 if (col_varExist(columnTrl,'RF')) then
2580 rainFlux(levelIndex,profileCount) = col_getElem(columnTrl,levelIndex,headerIndex,'RF')
2581 end if
2582 if (rainFlux(levelIndex,profileCount) < qlim_getMinValueCloud('RF') .or. &
2583 rainFlux(levelIndex,profileCount) > qlim_getMaxValueCloud('RF')) then
2584 write(*,*) 'tvs_fillProfiles: rainFlux=' , rainFlux(:,profileCount)
2585 call utl_abort('tvs_fillProfiles: columnTrl has rain flux outside RTTOV bounds')
2586 end if
2587
2588 ! snow flux (zero, if not part of control variables)
2589 if (col_varExist(columnTrl,'SF')) then
2590 snowFlux(levelIndex,profileCount) = col_getElem(columnTrl,levelIndex,headerIndex,'SF')
2591 end if
2592 if (snowFlux(levelIndex,profileCount) < qlim_getMinValueCloud('SF') .or. &
2593 snowFlux(levelIndex,profileCount) > qlim_getMaxValueCloud('SF')) then
2594 write(*,*) 'tvs_fillProfiles: snowFlux=' , snowFlux(:,profileCount)
2595 call utl_abort('tvs_fillProfiles: columnTrl has snow flux outside RTTOV bounds')
2596 end if
2597
2598 ! cloud fraction (zero, if not part of control variables)
2599 if (col_varExist(columnTrl,'CLDR')) then
2600 cloudFraction(levelIndex,profileCount) = col_getElem(columnTrl,levelIndex,headerIndex,'CLDR')
2601 end if
2602 if (cloudFraction(levelIndex,profileCount) < qlim_getMinValueCloud('CLDR') .or. &
2603 cloudFraction(levelIndex,profileCount) > qlim_getMaxValueCloud('CLDR')) then
2604 write(*,*) 'tvs_fillProfiles: cloudFraction=' , cloudFraction(:,profileCount)
2605 call utl_abort('tvs_fillProfiles: columnTrl has cloud fraction outside RTTOV bounds')
2606 end if
2607
2608 ciw(levelIndex,profileCount) = ciw(levelIndex,profileCount) * tvs_cloudScaleFactor
2609 rainFlux(levelIndex,profileCount) = rainFlux(levelIndex,profileCount) * tvs_cloudScaleFactor
2610 snowFlux(levelIndex,profileCount) = snowFlux(levelIndex,profileCount) * tvs_cloudScaleFactor
2611
2612 end if ! runObsOperatorWithHydrometeors .and. surfTypeIsWater
2613 end do ! levelIndex
2614
2615 if (tvs_coefs(sensorIndex) %coef % nozone > 0 .and. .not. tvs_useO3Climatology) then
2616 column_ptr => col_getColumn(columnTrl, headerIndex, trim(ozoneVarName) )
2617 ! Conversion from microgram/km to ppmv (to have the same units as climatology when tvs_useO3Climatology is .true.
2618 ! Conversion to kg/kg for use by RTTOV in done later
2619 ozone(:,profileCount) = column_ptr(:) * 1.0D-9 * o3Mixratio2ppmv
2620 end if
2621
2622 end do bobs2
2623
2624 ! 2.5 Get ozone profiles (ppmv) from climatology if necessary
2625 if (tvs_coefs(sensorIndex) %coef % nozone > 0 .and. tvs_useO3Climatology) then
2626 call ozo_get_profile (ozone, latitudes, pressure, nlv_T, profileCount)
2627 end if
2628
2629 ! 2.5 Fill profiles structure
2630
2631 do profileIndex = 1 , profileCount
2632 tovsIndex = sensorTovsIndexes(profileIndex)
2633 headerIndex = sensorHeaderIndexes(profileIndex)
2634 profiles(tovsIndex) % gas_units = gas_unit_specconc ! all gas profiles are supposed to be provided in kg/kg (specific humidity, i.e. mass mixing ratio [kg/kg] over wet air)
2635 profiles(tovsIndex) % id = '' ! profile id, up to 128 characters, to consider for use
2636 profiles(tovsIndex) % nlevels = nlv_T
2637 profiles(tovsIndex) % nlayers = nlv_T - 1
2638 profiles(tovsIndex) % date(1) = year
2639 profiles(tovsIndex) % date(2) = month
2640 profiles(tovsIndex) % date(3) = day
2641 profiles(tovsIndex) % latitude = latitudes(profileIndex)
2642 profiles(tovsIndex) % elevation = 0.001d0 * col_getHeight(columnTrl,ilowlvl_T,headerIndex,'TH') ! unite km
2643 call validateRttovVariable(profiles(tovsIndex) % elevation, "elevation", obsSpaceData, headerIndex, maximum_opt=elevmax)
2644 profiles(tovsIndex) % skin % watertype = watertype_ocean_water !utilise pour calcul rayonnement solaire reflechi seulement
2645 profiles(tovsIndex) % skin % t = col_getElem(columnTrl,1,headerIndex,'TG')
2646 call validateRttovVariable( profiles(tovsIndex) % skin % t, "skin temperature", obsSpaceData, headerIndex, tmin, tmax)
2647 profiles(tovsIndex) % skin % salinity = 35.d0 ! for FASTEM-4 only to revise (practical salinity units)
2648 profiles(tovsIndex) % skin % fastem(:) = 0.0d0
2649 profiles(tovsIndex) % skin % snow_fraction = 0.d0 ! Surface coverage snow fraction(0-1), used only by IR emissivity atlas
2650 profiles(tovsIndex) % skin % soil_moisture = 0.d0 ! soil moisure (m**3/m**3) not yet used
2651 profiles(tovsIndex) % s2m % t = col_getElem(columnTrl,ilowlvl_T,headerIndex,'TT')
2652 call validateRttovVariable(profiles(tovsIndex) % s2m % t, '2m air temperature', &
2653 obsSpaceData, headerIndex, tmin, tmax)
2654 profiles(tovsIndex) % s2m % p = col_getElem(columnTrl,1 ,headerIndex,'P0')*MPC_MBAR_PER_PA_R8
2655 call validateRttovVariable(profiles(tovsIndex) % s2m % p, 'surface pressure', &
2656 obsSpaceData, headerIndex, pmin, pmax)
2657 profiles(tovsIndex) % s2m % u = col_getElem(columnTrl,ilowlvl_M,headerIndex,'UU')
2658 profiles(tovsIndex) % s2m % v = col_getElem(columnTrl,ilowlvl_M,headerIndex,'VV')
2659 wind = sqrt( profiles(tovsIndex) % s2m % u ** 2 + &
2660 profiles(tovsIndex) % s2m % v ** 2 )
2661 if ( wind > wmax ) then
2662 write(*,*) 'tvs_fillProfiles: !!! WARNING !!!'
2663 write(*,*) 'tvs_fillProfiles: INVALID 10m wind speed'
2664 write(*,*) 'tvs_fillProfiles: headerIndex ', headerIndex, " !"
2665 write(*,*) 'tvs_fillProfiles: modulus ', wind, ' larger than ', wmax, 'set to zero !'
2666 profiles(tovsIndex) % s2m % u = 0.d0
2667 profiles(tovsIndex) % s2m % v = 0.d0
2668 call rejectObs(obsSpaceData, headerIndex)
2669 end if
2670 profiles(tovsIndex) % s2m % o = 0.0d0 !surface ozone never used
2671 profiles(tovsIndex) % s2m % wfetc = 100000.0d0 ! Wind fetch (in meter for rttov10 ?) used to calculate reflection of solar radiation by sea surface
2672 profiles(tovsIndex) % icede_param = 0
2673 profiles(tovsIndex) % Be = 0.4d0 ! earth magnetic field strength (gauss) (must be non zero)
2674 profiles(tovsIndex) % cosbk = 0.0d0 ! cosine of the angle between the earth magnetic field and wave propagation direction
2675 profiles(tovsIndex) % p(:) = pressure(:,profileIndex)
2676 call validateRttovVariable(profiles(tovsIndex) % p(nlv_T), "pressure profile near surface", obsSpaceData, headerIndex, maximum_opt=2000.d0)
2677 call validateRttovVariable(profiles(tovsIndex) % p(1), "pressure profile near top", obsSpaceData, headerIndex, 0.d0)
2678 !RTTOV scatt needs half pressure levels (see figure 5 of RTTOV 12 User's Guide)
2679 if (runObsOperatorWithHydrometeors) then
2680 cld_profiles(tovsIndex) % ph (1) = 0.d0
2681 cld_profiles(tovsIndex) % cfrac = 0.d0
2682 do levelIndex = 1, nlv_T - 1
2683 cld_profiles(tovsIndex) % ph (levelIndex+1) = 0.5d0 * (profiles(tovsIndex) % p(levelIndex) + &
2684 profiles(tovsIndex) % p(levelIndex+1))
2685 end do
2686 cld_profiles(tovsIndex) % ph (nlv_T+1) = profiles(tovsIndex) % s2m % p
2687 end if
2688 column_ptr => col_getColumn(columnTrl, headerIndex,'TT' )
2689 profiles(tovsIndex) % t(:) = column_ptr(:)
2690 call validateRttovProfile(profiles(tovsIndex) % t, 'temparature', tmin, tmax, obsSpaceData, headerIndex)
2691 if (tvs_coefs(sensorIndex) %coef %nozone > 0) then
2692 profiles(tovsIndex) % o3(:) = ozone(:,profileIndex) * o3ppmv2Mixratio ! Climatology output is ppmv over dry air ! because atmosphere is very dry where there is significant absorption by ozone)
2693 if (.not. tvs_useO3Climatology) then
2694 profiles(tovsIndex) % s2m % o = col_getElem(columnTrl,ilowlvl_T,headerIndex,trim(ozoneVarName)) * 1.0d-9 ! Assumes model ozone in ug/kg
2695 end if
2696 call validateRttovProfile(profiles(tovsIndex) % o3, 'ozone', o3min, o3max, obsSpaceData, headerIndex)
2697 end if
2698
2699 column_ptr => col_getColumn(columnTrl, headerIndex,'HU' )
2700 profiles(tovsIndex) % q(:) = column_ptr(:)
2701 call validateRttovProfile(profiles(tovsIndex) % q, 'water vapor', qmin, qmax, obsSpaceData, headerIndex)
2702 profiles(tovsIndex) % ctp = 1013.25d0
2703 profiles(tovsIndex) % cfraction = 0.d0
2704 if (runObsOperatorWithClw) then
2705 profiles(tovsIndex) % clw(:) = clw(:,profileIndex)
2706 else if (runObsOperatorWithHydrometeors) then
2707 cld_profiles(tovsIndex) % hydro(:,1) = rainFlux(:,profileIndex)
2708 cld_profiles(tovsIndex) % hydro(:,2) = snowFlux(:,profileIndex)
2709 cld_profiles(tovsIndex) % hydro(:,4) = clw(:,profileIndex)
2710 cld_profiles(tovsIndex) % hydro(:,5) = ciw(:,profileIndex)
2711
2712 do levelIndex = 1, nlv_T
2713 if (cld_profiles(tovsIndex) % hydro(levelIndex,1) > qlim_getMinValueCloud('RF') .or. &
2714 cld_profiles(tovsIndex) % hydro(levelIndex,2) > qlim_getMinValueCloud('SF') .or. &
2715 cld_profiles(tovsIndex) % hydro(levelIndex,4) > qlim_getMinValueCloud('LWCR') .or. &
2716 cld_profiles(tovsIndex) % hydro(levelIndex,5) > qlim_getMinValueCloud('IWCR')) then
2717
2718 ! set to overcast cloud, if CLDR not part of control variables
2719 if (col_varExist(columnTrl,'CLDR')) then
2720 cld_profiles(tovsIndex) % hydro_frac(levelIndex,1) = cloudFraction(levelIndex,profileIndex)
2721 else
2722 cld_profiles(tovsIndex) % hydro_frac(levelIndex,1) = 1.0d0
2723 end if
2724 else
2725 cld_profiles(tovsIndex) % hydro_frac(levelIndex,1) = qlim_getMinValueCloud('CLDR')
2726 end if
2727
2728 end do ! levelIndex
2729 end if ! runObsOperatorWithClw
2730 end do ! profileIndex
2731
2732 deallocate (pressure, stat = allocStatus(2))
2733 deallocate (ozone, stat = allocStatus(3))
2734 deallocate (latitudes, stat = allocStatus(4))
2735 deallocate (sensorHeaderIndexes, stat = allocStatus(5))
2736 deallocate (sensorTovsIndexes, stat = allocStatus(6))
2737 if (tvs_coefs(sensorIndex) %coef %nozone > 0 .and. .not.tvs_useO3Climatology) then
2738 deallocate (ozone, stat = allocStatus(7))
2739 end if
2740 if ( allocated(clw) ) then
2741 deallocate (clw, stat= allocStatus(8))
2742 end if
2743 if ( allocated(ciw) ) then
2744 deallocate (ciw, stat= allocStatus(9))
2745 deallocate (rainFlux, stat= allocStatus(10))
2746 deallocate (snowFlux, stat= allocStatus(11))
2747 deallocate (cloudFraction, stat= allocStatus(12))
2748 end if
2749 deallocate (surfTypeIsWater,stat= allocStatus(13))
2750
2751 call utl_checkAllocationStatus(allocStatus, ' tvs_fillProfiles', .false.)
2752
2753 end do sensor_loop
2754
2755 end subroutine tvs_fillProfiles
2756
2757 !--------------------------------------------------------------------------
2758 ! tvs_getCorrectedSatelliteAzimuth
2759 !--------------------------------------------------------------------------
2760 function tvs_getCorrectedSatelliteAzimuth(obsSpaceData, headerIndex) result(correctedAzimuth)
2761 !
2762 !:Purpose: get properly corrected satellite Azimuth Angle from obsSpaceData header
2763 !
2764 implicit none
2765
2766 ! Arguments:
2767 type(struct_obs), intent(in) :: obsSpaceData ! obsSpaceData structure
2768 integer, intent(in) :: headerIndex ! location in header
2769 ! Result:
2770 real(8) :: correctedAzimuth ! corrected azimuth (function result)
2771
2772 ! Locals:
2773 integer :: sensorNo, tovsIndex
2774
2775 correctedAzimuth = obs_headElem_r(obsSpaceData,OBS_AZA,headerIndex)
2776
2777 tovsIndex = tvs_tovsIndex (headerIndex)
2778 if ( tovsIndex < 0) return
2779
2780 sensorNo = tvs_lsensor(tovsIndex)
2781
2782 if ( .not. tvs_isAzimuthValid(sensorNo) ) then
2783 correctedAzimuth = obs_missingValue_R
2784 return
2785 end if
2786
2787 if ( tvs_doAzimuthCorrection(sensorNo) ) then
2788 ! Correction sur la definition de l'angle.
2789 correctedAzimuth = obs_headElem_r(obsSpaceData,OBS_SAZ,headerIndex) + correctedAzimuth
2790 if ( correctedAzimuth > 360.d0 ) correctedAzimuth = correctedAzimuth - 360.d0
2791 end if
2792
2793 end function tvs_getCorrectedSatelliteAzimuth
2794
2795
2796 !--------------------------------------------------------------------------
2797 ! tvs_rttov
2798 !--------------------------------------------------------------------------
2799 subroutine tvs_rttov(obsSpaceData, bgckMode, beSilent)
2800 !
2801 !:Purpose: Interface for RTTOV non linear operator
2802 ! tvs_fillProfiles should be called before
2803 !
2804 implicit none
2805
2806 ! Arguments:
2807 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData structure
2808 logical, intent(in) :: bgckMode ! flag to transfer transmittances and cloudy overcast radiances in bgck mode
2809 logical, intent(in) :: beSilent ! flag to control verbosity
2810
2811 ! Locals:
2812 integer :: nlv_T
2813 integer :: btCount
2814 integer :: allocStatus(4)
2815 integer :: rttov_err_stat ! rttov error return code
2816 integer, external :: omp_get_num_threads
2817 integer :: nthreads,max_nthreads
2818 integer :: sensorId, tovsIndex
2819 integer :: channelIndex, channelIndexFound, channelNumber
2820 integer :: profileCount
2821 integer :: profileIndex, levelIndex, jj, btIndex
2822 integer :: instrum
2823 integer :: sensorType ! sensor type (1=infrared; 2=microwave; 3=high resolution,4=polarimetric)
2824 integer, allocatable:: sensorTovsIndexes(:)
2825 type (rttov_emissivity), pointer :: emissivity_local(:) ! emissivity structure with input and output
2826 type (rttov_chanprof), pointer :: chanprof(:)
2827 type (rttov_chanprof), allocatable :: chanprof1(:)
2828 type (rttov_radiance) :: radiancedata_d, radiancedata_d1
2829 type (rttov_transmission) :: transmission, transmission1
2830 integer :: asw
2831 logical, pointer :: calcemis(:)
2832 real(8), allocatable :: surfem1(:)
2833 integer, allocatable :: frequencies(:)
2834 real(8), allocatable :: uOfWLandWSurfaceEmissivity(:)
2835 logical, allocatable :: lchannel_subset(:,:)
2836 integer :: profileIndex2, tb1, tb2
2837 integer :: istart, iend, bodyIndex, headerIndex
2838 real(8) :: clearMwRadiance
2839 logical :: ifBodyIndexFound
2840 logical :: runObsOperatorWithClw
2841 logical :: runObsOperatorWithHydrometeors
2842
2843 if ( .not. beSilent ) write(*,*) 'tvs_rttov: Starting'
2844 if ( .not. beSilent ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2845
2846 if (tvs_nobtov == 0) return ! exit if there are not tovs data
2847
2848 ! 1. Get number of threads available and allocate memory for some variables
2849 !$omp parallel
2850 max_nthreads = omp_get_num_threads()
2851 !$omp end parallel
2852 allocStatus(:) = 0
2853 allocate(sensorTovsIndexes(tvs_nobtov),stat=allocStatus(1))
2854 call utl_checkAllocationStatus(allocStatus(1:1), ' tvs_rttov sensorTovsIndexes')
2855
2856
2857 ! 1.1 Read surface information
2858 if ( bgckMode ) call EMIS_READ_CLIMATOLOGY
2859
2860 ! 2. Computation of hx for tovs data only
2861
2862 ! Loop over all sensors specified by user
2863 sensor_loop:do sensorId = 1, tvs_nsensors
2864
2865 sensorType = tvs_coefs(sensorId) % coef % id_sensor
2866 instrum = tvs_coefs(sensorId) % coef % id_inst
2867
2868 runObsOperatorWithClw = (tvs_numMWInstrumUsingCLW /= 0 .and. &
2869 tvs_isInstrumUsingCLW(tvs_instruments(sensorId)))
2870 runObsOperatorWithHydrometeors = (tvs_numMWInstrumUsingHydrometeors /= 0 .and. &
2871 tvs_isInstrumUsingHydrometeors(tvs_instruments(sensorId)))
2872
2873 if (runObsOperatorWithClw .and. runObsOperatorWithHydrometeors) then
2874 call utl_abort('tvs_rttov: this instrument is mentioned in using CLW and hydrometeors.')
2875 end if
2876
2877 ! loop over all obs.
2878 profileCount = 0
2879 obs_loop: do tovsIndex = 1, tvs_nobtov
2880
2881 ! Currently processed sensor?
2882 if ( tvs_lsensor(tovsIndex) == sensorId ) then
2883 profileCount = profileCount + 1
2884 sensorTovsIndexes(profileCount) = tovsIndex
2885 nlv_T = tvs_profiles_nl(tovsIndex) % nlevels
2886 end if
2887 end do obs_loop
2888
2889 if (profileCount == 0) cycle sensor_loop
2890
2891 ! 2.1 Calculate the actual number of threads which will be used.
2892
2893 nthreads = min(max_nthreads, profileCount )
2894
2895 ! 2.2 Prepare all input variables required by rttov.
2896
2897 if ( bgckMode .and. tvs_isInstrumHyperSpectral(instrum) ) then
2898 btCount = profileCount * tvs_nchan(sensorId)
2899 else
2900 btCount = tvs_countRadiances(sensorTovsIndexes(1:profileCount), obsSpaceData)
2901 end if
2902
2903 if (allocated(tvs_bodyIndexFromBtIndex)) deallocate(tvs_bodyIndexFromBtIndex)
2904 allocate(tvs_bodyIndexFromBtIndex(btCount))
2905
2906 if ( btCount == 0 ) cycle sensor_loop
2907 tvs_bodyIndexFromBtIndex(:) = -1
2908 allocate ( surfem1 (btCount) ,stat=allocStatus(1))
2909
2910 asw = 1 ! Allocation
2911 call rttov_alloc_direct( &
2912 allocStatus(2), &
2913 asw, &
2914 nprofiles=profileCount, & ! (not used)
2915 nchanprof=btCount, &
2916 nlevels=nlv_T, &
2917 chanprof=chanprof, &
2918 opts=tvs_opts(sensorId), &
2919 coefs=tvs_coefs(sensorId), &
2920 transmission=transmission, &
2921 radiance=radiancedata_d, &
2922 calcemis=calcemis, &
2923 emissivity=emissivity_local,&
2924 init=.true.)
2925
2926 if (runObsOperatorWithHydrometeors) then
2927 allocate ( frequencies(btCount), stat=allocStatus(3))
2928 end if
2929
2930 if (useUofWIREmiss) then
2931 allocate ( uOfWLandWSurfaceEmissivity(btCount), stat=allocStatus(4) )
2932 end if
2933 call utl_checkAllocationStatus(allocStatus, ' tvs_rttov')
2934
2935 ! get Hyperspectral IR emissivities
2936 if ( tvs_isInstrumHyperSpectral(instrum) ) then
2937 surfem1(:) = 0.
2938 if ( bgckMode ) then
2939 call emis_getIrEmissivity (surfem1,tvs_nchan(sensorId),sensorId,profileCount,btCount,sensorTovsIndexes)
2940 else
2941 call tvs_getHIREmissivities(sensorTovsIndexes(1:profileCount), obsSpaceData, surfem1)
2942 end if
2943 end if
2944
2945 if ( bgckMode .and. tvs_isInstrumHyperSpectral(instrum) ) then
2946 btIndex = 0
2947 do profileIndex = 1 , profileCount
2948 do channelIndex = 1,tvs_nchan(sensorId)
2949 btIndex = btIndex + 1
2950 chanprof(btIndex)%prof = profileIndex
2951 chanprof(btIndex)%chan = channelIndex
2952 end do
2953 end do
2954
2955 do profileIndex = 1, profileCount
2956 headerIndex = tvs_headerIndex(sensorTovsIndexes(profileIndex))
2957 if (headerIndex > 0) then
2958 istart = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2959 iend = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + istart - 1
2960 do bodyIndex = istart, iend
2961 call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
2962 channelNumber, channelIndex )
2963 if (channelIndex > 0) then
2964 tvs_bodyIndexFromBtIndex((profileIndex-1)*tvs_nchan(sensorId)+channelIndex) = bodyIndex
2965 else
2966 write(*,*) 'tvs_rttov: strange channel number',channelNumber
2967 end if
2968 end do
2969 end if
2970 end do
2971
2972 else
2973 allocate( lchannel_subset(profileCount,tvs_nchan(sensorId)) )
2974 call tvs_getChanprof(sensorTovsIndexes(1:profileCount), obsSpaceData, chanprof, lchannel_subset_opt = lchannel_subset, iptobs_cma_opt = tvs_bodyIndexFromBtIndex)
2975 if (runObsOperatorWithHydrometeors) then
2976 call rttov_scatt_setupindex ( &
2977 rttov_err_stat, &
2978 profileCount, & ! number of profiles
2979 tvs_nchan(sensorId), & ! number of channels
2980 tvs_coefs(sensorId), & ! coef structure read in from rttov coef file
2981 tvs_coef_scatt(sensorId), & !
2982 btcount, & ! number of calculated channels
2983 chanprof, & ! channels and profile numbers
2984 frequencies, & ! array, frequency number for each channel
2985 lchannel_subset ) ! OPTIONAL array of logical flags to indicate a subset of channels
2986 end if
2987 deallocate( lchannel_subset )
2988 end if
2989
2990 call tvs_getOtherEmissivities(chanprof, sensorTovsIndexes, sensorType, instrum, surfem1, calcemis)
2991
2992 if (useUofWIREmiss .and. tvs_isInstrumHyperSpectral(instrum) .and. bgckMode) then
2993 if (.not. allocated (tvs_atlas)) allocate(tvs_atlas(tvs_nsensors))
2994 if ( .not. tvs_atlas(sensorId)%init) then
2995 call rttov_setup_emis_atlas( rttov_err_stat, &! out
2996 tvs_opts(sensorId), &! in
2997 tvs_profiles_nl(1)%date(2) , &! in
2998 atlas_type_ir, &! in
2999 tvs_atlas(sensorId), &! in
3000 ir_atlas_ang_corr = .false., &! in
3001 ir_atlas_read_std = .false., &! in
3002 coefs = tvs_coefs(sensorId) )
3003 if (rttov_err_stat/=0) then
3004 write(*,*) 'Error in rttov_atlas_setup IR',rttov_err_stat
3005 call utl_abort('tvs_rttov')
3006 end if
3007 end if
3008
3009 call rttov_get_emis( rttov_err_stat, & ! out
3010 tvs_opts(sensorId), & ! in
3011 chanprof(1:btCount), & ! in
3012 tvs_profiles_nl(sensorTovsIndexes(1:profileCount)), & ! in
3013 tvs_coefs(sensorId), & ! in
3014 tvs_atlas(sensorId), & ! inout
3015 uOfWLandWSurfaceEmissivity(1:btCount) ) ! out
3016
3017 if (rttov_err_stat /= 0) then
3018 write(*,*) 'Error in rttov_get_emis IR', rttov_err_stat
3019 call utl_abort('tvs_rttov')
3020 end if
3021
3022 do profileIndex=1, profileCount !loop on profiles
3023 jj = sensorTovsIndexes(profileIndex)
3024 do btIndex=1, btCount !loop on channels
3025 if (chanprof(btIndex)%prof==profileIndex) then
3026 ! surftype: 0 land, 1 sea, 2 sea-ice
3027 ! this logic is primitive and could be improved for example using
3028 ! additional criteria based on emissivity_std and emissivity_flg
3029 !Definition of emis_flag:
3030 ! emis_flag:Flag_0 = '0 = sea, no MOD11 data' ;
3031 ! emis_flag:Flag_1 = '1 = land where BF method was applied' ;
3032 ! emis_flag:Flag_2 = '2 = land where data was filled with average (original UWiremis bfemis_flag=2 or 3 or 4' ;
3033 ! emis_flag:Flag_3 = '3 = contains inland water or coastline by the sea/land mask where the BF method was used' ;
3034 ! emis_flag:Flag_4 = '4 = contains inland water or coastline by the sea/land mask where data was filled with average original UWiremis bfemis_flag=2 or 3 or 4' ;
3035 ! emis_flag:Flag_5 = '5 = contains coastline by land fraction where the BF method was used' ;
3036 ! emis_flag:Flag_6 = '6 = contains coastline by land fraction where data was filled with average (original UWiremis bfemis_flag=2 or 3 or 4' ;
3037 ! other information that could be useful for quality control can be found in the in the profile_qc structure
3038 ! Now we have the 'traditionnal' emissivity in surfem1(:)
3039 ! and University of Wisconsin emissivity in uOfWLandWSurfaceEmissivity(:)
3040 if (tvs_profiles_nl(jj)% skin % surftype == surftype_land .and. &
3041 uOfWLandWSurfaceEmissivity(btIndex) > 0.5 ) then
3042 emissivity_local(btIndex)%emis_in = uOfWLandWSurfaceEmissivity(btIndex)
3043 else
3044 emissivity_local(btIndex)%emis_in = surfem1(btIndex)
3045 end if
3046 end if
3047 end do
3048 end do
3049
3050 else if (sensorType == sensor_id_mw) then
3051 call tvs_getMWemissivityFromAtlas(surfem1(1:btcount), emissivity_local, sensorId, chanprof, &
3052 sensorTovsIndexes(1:profileCount))
3053 else
3054 emissivity_local(:)%emis_in = surfem1(:)
3055 end if
3056
3057 ! 2.3 Compute radiance with rttov_direct
3058
3059 rttov_err_stat = 0
3060
3061 if( bgckMode .and. tvs_isInstrumHyperSpectral(instrum) ) then
3062 write(*,*) 'for bgck IR: call rttov_parallel_direct for each profile...'
3063
3064 asw = 1 ! 1 to allocate,0 to deallocate
3065 ! allocate transmitance structure for 1 profile
3066 call rttov_alloc_transmission(allocStatus(1), transmission1, nlevels=nlv_T, &
3067 nchanprof=tvs_nchan(sensorId), asw=asw, init=.true.)
3068 ! allocate radiance structure for 1 profile
3069 call rttov_alloc_rad (allocStatus(2),tvs_nchan(sensorId), radiancedata_d1,nlv_T,asw,init=.true.)
3070 ! allocate chanprof for 1 profile
3071 allocate(chanprof1(tvs_nchan(sensorId)))
3072 do channelIndex = 1,tvs_nchan(sensorId)
3073 chanprof1(channelIndex)%prof = 1
3074 chanprof1(channelIndex)%chan = channelIndex
3075 end do
3076
3077 do profileIndex2 = 1, profileCount
3078 tb1 = 1 + (profileIndex2-1) * tvs_nchan(sensorId)
3079 tb2 = profileIndex2 * tvs_nchan(sensorId)
3080
3081 call rttov_parallel_direct( &
3082 rttov_err_stat, & ! out
3083 chanprof1, & ! in
3084 tvs_opts(sensorId), & ! in
3085 tvs_profiles_nl(sensorTovsIndexes(profileIndex2):sensorTovsIndexes(profileIndex2)), & ! in
3086 tvs_coefs(sensorId), & ! in
3087 transmission1, & ! inout
3088 radiancedata_d1, & ! inout
3089 calcemis=calcemis(tb1:tb2), & ! in
3090 emissivity=emissivity_local(tb1:tb2), & ! inout
3091 nthreads=nthreads )
3092
3093 ! copy contents of single profile structures into complete structures
3094 transmission%tau_total(tb1:tb2) = transmission1%tau_total(:)
3095 transmission%tau_levels(:,tb1:tb2) = transmission1%tau_levels(:,:)
3096 transmission%tausun_levels_path1(:,tb1:tb2) = transmission1%tausun_levels_path1(:,:)
3097 transmission%tausun_levels_path2(:,tb1:tb2) = transmission1%tausun_levels_path2(:,:)
3098 transmission%tausun_total_path1(tb1:tb2) = transmission1%tausun_total_path1(:)
3099 transmission%tausun_total_path2(tb1:tb2) = transmission1%tausun_total_path2(:)
3100 radiancedata_d%clear(tb1:tb2) = radiancedata_d1%clear(:)
3101 radiancedata_d%total(tb1:tb2) = radiancedata_d1%total(:)
3102 radiancedata_d%bt_clear(tb1:tb2) = radiancedata_d1%bt_clear(:)
3103 radiancedata_d%bt(tb1:tb2) = radiancedata_d1%bt(:)
3104 radiancedata_d%refl_clear(tb1:tb2) = radiancedata_d1%refl_clear(:)
3105 radiancedata_d%refl(tb1:tb2) = radiancedata_d1%refl(:)
3106 radiancedata_d%overcast(:,tb1:tb2) = radiancedata_d1%overcast(:,:)
3107 radiancedata_d%cloudy(tb1:tb2) = radiancedata_d1%cloudy(:)
3108
3109 end do
3110
3111 ! transmittance deallocation for 1 profile
3112 deallocate(chanprof1)
3113 asw = 0 ! 1 to allocate,0 to deallocate
3114 ! transmittance deallocation for 1 profile
3115 call rttov_alloc_transmission(allocStatus(1),transmission1,nlevels=nlv_T, &
3116 nchanprof=tvs_nchan(sensorId), asw=asw )
3117 ! radiance deallocation for 1 profile
3118 call rttov_alloc_rad (allocStatus(2), tvs_nchan(sensorId), radiancedata_d1, nlv_T, asw)
3119
3120 else
3121
3122 ! run clear-sky RTTOV, save the radiances in OBS_BTCL of obsSpaceData
3123 if ((runObsOperatorWithClw .or. runObsOperatorWithHydrometeors) .and. &
3124 obs_columnActive_RB(obsSpaceData, OBS_BTCL)) then
3125
3126 ! run rttovScatt
3127 if (runObsOperatorWithHydrometeors) then
3128 ! set the cloud profile in tvs_cld_profiles_nl to zero
3129 call updateCloudInTovsCloudProfile(sensorTovsIndexes(1:profileCount), &
3130 nlv_T, &
3131 mode='save', &
3132 beSilent=.true.)
3133 call rttov_scatt( &
3134 rttov_err_stat, &! out
3135 tvs_opts_scatt(sensorId), &! in
3136 nlv_T, &! in
3137 chanprof, &! in
3138 frequencies, &! in
3139 tvs_profiles_nl(sensorTovsIndexes(1:profileCount)), &! in
3140 tvs_cld_profiles_nl(sensorTovsIndexes(1:profileCount)), &! in
3141 tvs_coefs(sensorId), &! in
3142 tvs_coef_scatt(sensorId), &! in
3143 calcemis, &! in
3144 emissivity_local, &! inout
3145 radiancedata_d)
3146 else
3147 ! set the cloud profile in tvs_profiles_nl to zero
3148 call updateCloudInTovsProfile(sensorTovsIndexes(1:profileCount), &
3149 nlv_T, &
3150 mode='save', &
3151 beSilent=.true.)
3152 call rttov_parallel_direct( &
3153 rttov_err_stat, & ! out
3154 chanprof, & ! in
3155 tvs_opts(sensorId), & ! in
3156 tvs_profiles_nl(sensorTovsIndexes(1:profileCount)), & ! in
3157 tvs_coefs(sensorId), & ! in
3158 transmission, & ! inout
3159 radiancedata_d, & ! inout
3160 calcemis=calcemis, & ! in
3161 emissivity=emissivity_local, & ! inout
3162 nthreads=nthreads )
3163 end if ! run rttovScatt
3164
3165 ! save in obsSpaceData
3166 loopClearSky1: do btIndex = 1, btCount
3167 profileIndex = chanprof(btIndex)%prof
3168 channelIndex = chanprof(btIndex)%chan
3169 tovsIndex = sensorTovsIndexes(profileIndex)
3170
3171 clearMwRadiance = radiancedata_d % bt(btIndex)
3172
3173 headerIndex = tvs_headerIndex(tovsIndex)
3174 if ( headerIndex < 1 ) cycle loopClearSky1
3175 istart = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
3176 iend = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) + istart - 1
3177
3178 ifBodyIndexFound = .false.
3179 loopClearSky2: do bodyIndex = istart, iend
3180 call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
3181 channelNumber, channelIndexFound )
3182 if ( channelIndex == channelIndexFound ) then
3183 ifBodyIndexFound = .true.
3184 exit loopClearSky2
3185 end if
3186 end do loopClearSky2
3187
3188 if ( .not. ifBodyIndexFound ) call utl_abort('tvs_rttov: bodyIndex not found.')
3189
3190 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
3191 call obs_bodySet_r(obsSpaceData, OBS_BTCL, bodyIndex, clearMwRadiance)
3192 end if
3193 end do loopClearSky1
3194
3195 ! restore the cloud profiles in ...
3196 if (runObsOperatorWithHydrometeors) then
3197 ! tvs_cld_profiles_nl
3198 call updateCloudInTovsCloudProfile(sensorTovsIndexes(1:profileCount), &
3199 nlv_T, &
3200 mode='restore', &
3201 beSilent=.true.)
3202 else
3203 ! tvs_profiles_nl
3204 call updateCloudInTovsProfile(sensorTovsIndexes(1:profileCount), &
3205 nlv_T, &
3206 mode='restore', &
3207 beSilent=.true.)
3208 end if
3209 end if ! run clear-sky RTTOV
3210
3211 if (runObsOperatorWithHydrometeors) then
3212 if (.not. beSilent) write(*,*) 'before rttov_scatt...', sensorID, profileCount
3213 call rttov_scatt( &
3214 rttov_err_stat, &! out
3215 tvs_opts_scatt(sensorId), &! in
3216 nlv_T, &! in
3217 chanprof, &! in
3218 frequencies, &! in
3219 tvs_profiles_nl(sensorTovsIndexes(1:profileCount)), &! in
3220 tvs_cld_profiles_nl(sensorTovsIndexes(1:profileCount)), &! in
3221 tvs_coefs(sensorId), &! in
3222 tvs_coef_scatt(sensorId), &! in
3223 calcemis, &! in
3224 emissivity_local, &! inout
3225 radiancedata_d)
3226 if ( .not. beSilent ) write(*,*) 'after rttov_scatt...'
3227 else
3228 if (.not. beSilent) write(*,*) 'before rttov_parallel_direct...', sensorID, profileCount
3229
3230 call rttov_parallel_direct( &
3231 rttov_err_stat, & ! out
3232 chanprof, & ! in
3233 tvs_opts(sensorId), & ! in
3234 tvs_profiles_nl(sensorTovsIndexes(1:profileCount)), & ! in
3235 tvs_coefs(sensorId), & ! in
3236 transmission, & ! inout
3237 radiancedata_d, & ! inout
3238 calcemis=calcemis, & ! in
3239 emissivity=emissivity_local, & ! inout
3240 nthreads=nthreads )
3241 if ( .not. beSilent ) write(*,*) 'after rttov_parallel_direct...'
3242
3243 end if
3244
3245 end if ! if (bgckMode .and. tvs_isInstrumHyperSpectral(instrum))
3246
3247 if ( .not. beSilent ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
3248
3249 if (rttov_err_stat /= 0) then
3250 write(*,*) 'Error in rttov_parallel_direct',rttov_err_stat
3251 call utl_abort('tvs_rttov')
3252 end if
3253
3254 ! 2.4 Store hx in the structure tvs_radiance
3255
3256 do btIndex = 1, btCount
3257 profileIndex = chanprof(btIndex)%prof
3258 channelIndex = chanprof(btIndex)%chan
3259 tovsIndex = sensorTovsIndexes(profileIndex)
3260 tvs_radiance(tovsIndex) % bt(channelIndex) = radiancedata_d % bt(btIndex)
3261
3262 if ( bgckMode ) then
3263 if ( .not. associated(tvs_radiance(tovsIndex) % clear)) then
3264 allocStatus = 0
3265 allocate( tvs_radiance(tovsIndex) % clear ( tvs_nchan(sensorId) ), stat= allocStatus(1) )
3266 ! allocate overcast black cloud sky radiance output
3267 allocate( tvs_radiance(tovsIndex) % overcast (nlv_T - 1, tvs_nchan(sensorId) ), stat=allocStatus(2))
3268 call utl_checkAllocationStatus(allocStatus(1:2), ' tvs_rttov')
3269 end if
3270 tvs_radiance(tovsIndex) % clear(channelIndex) = &
3271 radiancedata_d %clear(btIndex)
3272 do levelIndex = 1, nlv_T - 1
3273 tvs_radiance(tovsIndex) % overcast(levelIndex,channelIndex) = &
3274 radiancedata_d % overcast(levelIndex,btIndex)
3275 end do
3276 if (.not. allocated(tvs_transmission)) call tvs_allocTransmission(nlv_T)
3277 end if
3278
3279 if ( allocated( tvs_transmission) ) then
3280 do levelIndex = 1, nlv_T
3281 tvs_transmission(tovsIndex) % tau_levels(levelIndex,channelIndex) = &
3282 transmission % tau_levels(levelIndex,btIndex)
3283 end do
3284
3285 tvs_transmission(tovsIndex) % tau_total(channelIndex) = &
3286 transmission % tau_total(btIndex)
3287 end if
3288
3289 if ( allocated(tvs_emissivity) ) then
3290 tvs_emissivity(channelIndex,tovsIndex) = emissivity_local(btIndex)%emis_out
3291 end if
3292
3293 end do
3294
3295 ! Append Surface Emissivity into ObsSpaceData
3296 do btIndex = 1, btCount
3297 bodyIndex = tvs_bodyIndexFromBtIndex(btIndex)
3298 if (bodyIndex > 0) call obs_bodySet_r(obsSpaceData, OBS_SEM, bodyIndex, emissivity_local(btIndex)%emis_out)
3299 end do
3300
3301 ! Deallocate memory
3302 asw = 0 ! 0 to deallocate
3303 call rttov_alloc_direct( &
3304 allocStatus(1), &
3305 asw, &
3306 nprofiles=profileCount, & ! (not used)
3307 nchanprof=btCount, &
3308 nlevels=nlv_T, &
3309 chanprof=chanprof, &
3310 opts=tvs_opts(sensorId), &
3311 coefs=tvs_coefs(sensorId), &
3312 transmission=transmission, &
3313 radiance=radiancedata_d, &
3314 calcemis=calcemis, &
3315 emissivity=emissivity_local,&
3316 init=.true.)
3317
3318
3319 if (useUofWIREmiss) then
3320 deallocate ( uOfWLandWSurfaceEmissivity ,stat=allocStatus(2) )
3321 end if
3322 deallocate ( surfem1 ,stat=allocStatus(3) )
3323 if (allocated(frequencies)) deallocate (frequencies, stat=allocStatus(4))
3324 call utl_checkAllocationStatus(allocStatus, ' tvs_rttov', .false.)
3325
3326 end do sensor_loop
3327
3328 deallocate(sensorTovsIndexes)
3329
3330 end subroutine tvs_rttov
3331
3332 !--------------------------------------------------------------------------
3333 ! tvs_getMWemissivityFromAtlas
3334 !--------------------------------------------------------------------------
3335 subroutine tvs_getMWemissivityFromAtlas(originalEmissivity, updatedEmissivity, sensorId, chanprof, sensorTovsIndexes)
3336 implicit none
3337
3338 ! Arguments:
3339 real(8), intent(in) :: originalEmissivity(:)
3340 type(rttov_emissivity), intent(out) :: updatedEmissivity(:)
3341 integer, intent(in) :: sensorId
3342 type(rttov_chanprof), intent(in) :: chanprof(:)
3343 integer, intent(in) :: sensorTovsIndexes(:)
3344
3345 ! Locals:
3346 integer :: returnCode
3347 real(8) :: mWAtlasSurfaceEmissivity(size(originalEmissivity))
3348 integer :: btCount, profileCount
3349 integer :: profileIndex, btIndex, sensorIndex
3350
3351 btCount = size( originalEmissivity )
3352 if (useMWEmissivityAtlas) then
3353
3354 if (.not. allocated (tvs_atlas)) allocate(tvs_atlas(tvs_nsensors))
3355 if ( .not. tvs_atlas(sensorId)%init) then
3356 call rttov_setup_emis_atlas( returnCode, &! out
3357 tvs_opts(sensorId), &! in
3358 tvs_profiles_nl(1)%date(2), &! in
3359 atlas_type_mw, &! in
3360 tvs_atlas(sensorId), &! inout
3361 atlas_id = mWAtlasId, &! in ! 1 TELSEM2, 2 CNRM
3362 coefs = tvs_coefs(sensorId) )
3363 if (returnCode /= 0) then
3364 write(*,*) 'Error in rttov_atlas_setup MW',returnCode
3365 call utl_abort('tvs_getMWemissivityFromAtlas')
3366 end if
3367 end if
3368
3369 call rttov_get_emis( returnCode, & ! out
3370 tvs_opts(sensorId), & ! in
3371 chanprof, & ! in
3372 tvs_profiles_nl(sensorTovsIndexes(:)), & ! in
3373 tvs_coefs(sensorId), & ! in
3374 tvs_atlas(sensorId), & ! in
3375 mWAtlasSurfaceEmissivity) ! out
3376
3377 if (returnCode /= 0) then
3378 write(*,*) 'Error in rttov_get_emis MW', returnCode
3379 call utl_abort('tvs_getMWemissivityFromAtlas')
3380 end if
3381
3382 profileCount = size( sensorTovsIndexes )
3383
3384 do profileIndex=1, profileCount !loop on profiles
3385 sensorIndex = sensorTovsIndexes(profileIndex)
3386
3387 do btIndex=1, btCount !loop on channels
3388 if (chanprof(btIndex)%prof==profileIndex) then
3389 ! Now we have 0.75 in originalEmissivity(:) for land and sea ice
3390 ! and the MW atlas emissivity in mWAtlasSurfaceEmissivity(:)
3391
3392 if ( tvs_profiles_nl(sensorIndex)% skin % surftype == surftype_land .and. &
3393 mWAtlasSurfaceEmissivity(btIndex) > 0.d0 .and. &
3394 mWAtlasSurfaceEmissivity(btIndex) <= 1.d0 ) then ! check for missing values
3395 updatedEmissivity(btIndex)%emis_in = mWAtlasSurfaceEmissivity(btIndex)
3396
3397 else
3398 updatedEmissivity(btIndex)%emis_in = originalEmissivity(btIndex)
3399 end if
3400 ! Note that emissivity above sea-ice is not modified
3401 end if
3402 end do
3403 end do
3404 else
3405 updatedEmissivity(:)%emis_in = originalEmissivity(:)
3406 end if
3407 end subroutine tvs_getMWemissivityFromAtlas
3408
3409 !--------------------------------------------------------------------------
3410 ! comp_ir_emiss
3411 !--------------------------------------------------------------------------
3412 subroutine comp_ir_emiss (emiss, wind, angle, nchn, np, mchannel)
3413 !
3414 !:Purpose: Computes water infrared emissivity for a specific set of
3415 ! channel indices, wind speed and zenith angle.
3416 !
3417 implicit none
3418
3419 ! Arguments:
3420 real(8), intent(out) :: emiss(nchn,np) ! emissivities (0.-1.)
3421 real(8), intent(in) :: wind(np) ! wind: surface wind speed (m/s)
3422 real(8), intent(in) :: angle(np) ! angle: viewing angle (deg)
3423 integer, intent(in) :: nchn ! number of channels to process
3424 integer, intent(in) :: np !number of locations
3425 integer, intent(in) :: mchannel(nchn) ! vector of channel indices to process
3426
3427 ! Locals:
3428 integer, parameter :: MaxWn = 19
3429 integer, parameter :: Nparm=3
3430 integer, parameter :: MaxChan=19
3431
3432 real (8),parameter :: Theta(Nparm,MaxWn) = (/ &
3433 1700.381d0, 25.28534d0, 144.1023d0, &
3434 1738.149d0, 25.67787d0, 146.6139d0, &
3435 1769.553d0, 26.05250d0, 148.6586d0, &
3436 1778.610d0, 26.12333d0, 149.5127d0, &
3437 1794.245d0, 26.18523d0, 150.5874d0, &
3438 1791.904d0, 26.19991d0, 150.7076d0, &
3439 1806.872d0, 26.37132d0, 151.7191d0, &
3440 1926.078d0, 27.63825d0, 160.7103d0, &
3441 1969.155d0, 28.02767d0, 163.6069d0, &
3442 1975.549d0, 27.86465d0, 164.6228d0, &
3443 1991.288d0, 27.94312d0, 166.2924d0, &
3444 2082.691d0, 28.93558d0, 172.4025d0, &
3445 2182.872d0, 29.89974d0, 179.5839d0, &
3446 2338.510d0, 31.27507d0, 191.2063d0, &
3447 2164.615d0, 28.97152d0, 182.6279d0, &
3448 2099.714d0, 29.91868d0, 178.4015d0, &
3449 1857.644d0, 29.13640d0, 160.9822d0, &
3450 1610.696d0, 26.48602d0, 142.2768d0, &
3451 1503.969d0, 24.97931d0, 133.4392d0 /)
3452
3453 real (8),parameter :: C(Nparm,2,MaxWn) = (/ &
3454 0.9715104043561414d0,-1.2034233230944147D-06, -5.8742655960993913D-07, &
3455 0.9263932848727608d0,-9.4908630939690859D-04, 2.2831134823358876D-05, &
3456 0.9732503924722753d0,-1.2007007329295099D-06, -5.8767355551283423D-07, &
3457 0.9290947860585505d0,-9.5233413988900161D-04, 2.2640835623043761D-05, &
3458 0.9745005204317289d0, 1.2857517639804244D-06, -7.1356127087301190D-07, &
3459 0.9310852809117095d0,-9.5453509182819095D-04, 2.2562638663187251D-05, &
3460 0.9756204829761132d0, 1.2979181109743976D-06, -7.1406809362820139D-07, &
3461 0.9329073568177888d0,-9.5627536945214183D-04, 2.2442358508999558D-05, &
3462 0.9764012672766408d0,-2.0826654381361387D-06, -4.9103920569405721D-07, &
3463 0.9341937281933334d0,-9.5764423928102976D-04, 2.2326701573603621D-05, &
3464 0.9770513558720460d0, 4.1867599593267133D-07, -6.1768073971231931D-07, &
3465 0.9352981872014672d0,-9.5833614545300181D-04, 2.2261996883974513D-05, &
3466 0.9775970810179080d0,-1.2289690625562906D-06, -5.2953762169985775D-07, &
3467 0.9362188153954743d0,-9.5950872922696905D-04, 2.2251301675423482D-05, &
3468 0.9830610391451819d0, 2.7693589475690676D-07, -5.1580217018207558D-07, &
3469 0.9461121192685766d0,-9.5718115604053031D-04, 2.1087308573177295D-05, &
3470 0.9840097930773377d0,-1.4987900189155091D-06, -3.8281408128977038D-07, &
3471 0.9479758694804105d0,-9.5451252460440695D-04, 2.0800627740862229D-05, &
3472 0.9851056150728580d0,-6.5768237152417477D-07, -4.2053769829400935D-07, &
3473 0.9502084544618980d0,-9.4965534997704157D-04, 2.0326602209199427D-05, &
3474 0.9862706396188835d0,-2.3713068057993353D-06, -2.8671134918457728D-07, &
3475 0.9526580467595886d0,-9.4614505430749598D-04, 2.0001856872595840D-05, &
3476 0.9875307221489201d0, 1.3003462826947714D-07, -4.1335288320283954D-07, &
3477 0.9554195617948236d0,-9.3806678196435643D-04, 1.9407754748128057D-05, &
3478 0.9891153260567763d0,-8.0730206675976713D-07, -3.1811320626834656D-07, &
3479 0.9590558393678170d0,-9.2716287670223167D-04, 1.8690586764925213D-05, &
3480 0.9913304557147524d0,-2.1153391229093421D-08, -3.1094269595901165D-07, &
3481 0.9644162604969492d0,-9.0342753739935612D-04, 1.7274993357160937D-05, &
3482 0.9925188366950193d0,-4.6365959315123653D-07, -2.7020120347068712D-07, &
3483 0.9667877170960085d0,-9.0665804037922043D-04, 1.7083616616646458D-05, &
3484 0.9919408379810360d0,-2.0563508815953840D-06, -1.8066722718403761D-07, &
3485 0.9627535343397309d0,-9.7537134133678965D-04, 1.9698263973541952D-05, &
3486 0.9889406296815972d0,-2.3713068057993353D-06, -2.8671134918457728D-07, &
3487 0.9506051906192242d0,-1.0642902225813857D-03, 2.4235485973033298D-05, &
3488 0.9828819693848310d0,-7.4086701870172759D-07, -6.2949258820534062D-07, &
3489 0.9329616683158125d0,-1.0728027288012200D-03, 2.7209071863380586D-05, &
3490 0.9767410313266288d0,-9.1750038410238915D-07, -7.9177921107781349D-07, &
3491 0.9192775350344998d0,-1.0369254272157462D-03, 2.8000594542037504D-05 /)
3492
3493 real (8) a(MaxChan),b(MaxChan),cc(MaxChan) ! local variable
3494 real (8) ww
3495 integer Index,Ichan,IP
3496
3497 do Ichan = 1 , Nchn
3498
3499 Index = Mchannel(Ichan)
3500
3501 do IP=1,NP
3502
3503 ww = wind(IP)
3504 a(Ichan) = c(1,1,Index) + c(2,1,Index) * ww &
3505 + c(3,1,Index) * ww * ww
3506 b(Ichan) = c(1,2,Index) + c(2,2,Index) * ww &
3507 + c(3,2,Index)* ww * ww
3508
3509 cc(Ichan) = Theta(1,Index) + Theta(2,Index) * ww
3510
3511 emiss(Ichan,IP) = a(Ichan) + (b(Ichan) - a(Ichan)) * &
3512 exp(( (theta(3,Index) - 60.d0)**2.d0 &
3513 - (angle(IP) - theta(3,Index))**2.d0 ) / CC(Ichan))
3514
3515 end do
3516
3517 end do
3518
3519 end subroutine comp_ir_emiss
3520
3521 !--------------------------------------------------------------------------
3522 ! pcnt_box
3523 !--------------------------------------------------------------------------
3524 subroutine pcnt_box(f_low, f_high, nprf, ilat, ilon, klat, klon, ireduc)
3525 !
3526 !:Purpose: Computes a low resolution feature form a high
3527 ! resolution one by averaging.
3528 ! example: use for percentage of water
3529 !
3530 implicit none
3531
3532 ! Arguments:
3533 real(8), intent(out) :: f_low(nprf) ! Low resolution field
3534 real(8), intent(in) :: f_high(klon, klat)! High resolution field
3535 integer, intent(in) :: nprf ! Number of profiles
3536 integer, intent(in) :: ilat(nprf) ! Y-coordinate of profile
3537 integer, intent(in) :: ilon(nprf) ! X-coordinate of profile
3538 integer, intent(in) :: klon ! Max value of latitude indices
3539 integer, intent(in) :: klat ! Max value of longitude indices
3540 integer, intent(in) :: ireduc ! Means a 2xireduc+1 by 2xireduc+1 averaging
3541
3542 ! Locals
3543 integer :: nplon, jdlo1, jdlo2, jlon1, jlon2
3544 integer :: nx, ilat1, ilat2, ilon1, ilon2, jn, ii, jj
3545
3546 profiles : do jn = 1,nprf
3547
3548 nplon = 0
3549
3550 ! normal limits
3551
3552 ilat1=max(ilat(jn)-ireduc,1)
3553 ilat2=min(ilat(jn)+ireduc,klat)
3554 ilon1=max(ilon(jn)-ireduc,1)
3555 ilon2=min(ilon(jn)+ireduc,klon)
3556
3557 if (ilon1 == 1 .or. ilon2 == klon) then
3558 ! border cases for longitudes
3559 jdlo1 = ilon(jn)-ireduc
3560 jdlo2 = ilon(jn)+ireduc
3561
3562 if ( jdlo1 <= 0 ) then
3563 nplon = 1
3564 jlon1 = klon + jdlo1
3565 jlon2 = klon
3566 else if ( jdlo2 > klon ) then
3567 nplon = 1
3568 jlon1 = 1
3569 jlon2 = jdlo2 - klon
3570 end if
3571 end if
3572
3573 nx = 0
3574 f_low(jn) = 0.d0
3575
3576 do jj = ilat1, ilat2
3577
3578 do ii = ilon1, ilon2
3579 nx = nx + 1
3580 f_low(jn) = f_low(jn) + f_high(ii,jj)
3581 end do
3582
3583 if (nplon == 1) then
3584 ! additional cases at border 1-klon
3585 do ii = jlon1, jlon2
3586 nx = nx + 1
3587 f_low(jn) = f_low(jn) + f_high(ii,jj)
3588 end do
3589 end if
3590
3591 end do
3592
3593 f_low(jn) = f_low(jn) / dble(nx)
3594
3595 end do profiles
3596
3597 end subroutine pcnt_box
3598
3599 !--------------------------------------------------------------------------
3600 ! emis_read_climatology
3601 !--------------------------------------------------------------------------
3602 subroutine emis_read_climatology
3603 !
3604 !:Purpose: Read information about ceres surface type and water fraction.
3605 !
3606 implicit none
3607
3608 ! Locals:
3609 integer :: nisf,njsf,nksf
3610 integer :: niwa,njwa,nkwa
3611 character(len=20) :: ceresFile
3612 integer, external :: fnom,fstouv,fstfrm,fclos,fstlir
3613 integer :: isftest
3614 integer :: iv1,iv2,iv3,iv4,iv5,iv6
3615
3616 isftest = 0
3617
3618 !* get surface type and water fraction
3619 ceresFile = 'ceres_global.std'
3620 iv1 = fnom(isftest,ceresFile,'RND+R/O',0)
3621 iv2 = fstouv(isftest,'RND')
3622 iv3 = fstlir(surfaceType,isftest,nisf,njsf,nksf,-1,'SFC-TYPE',-1,-1,-1,'','TY')
3623 iv4 = utl_fstlir(waterFraction,isftest,niwa,njwa,nkwa,-1,'WATER_FR',-1,-1,-1,'','W%')
3624 iv5 = fstfrm(isftest)
3625 iv6 = fclos(isftest)
3626
3627 if (iv1 < 0 .or. iv2 < 0 .or. iv3 < 0 .or. iv4 < 0 .or. iv5 < 0 .or. iv6 < 0) then
3628 write(*,*) 'LES iv DE CERES ',iv1,iv2,iv3,iv4,iv5,iv6
3629 write(*,*) 'THESE NUMBER SHOULD NOT BE NEGATIVE WHEN DOING AIRS BACKGROUND CHECK'
3630 call utl_abort('Problem with file ceres_global.std in emis_read_climatology ')
3631 end if
3632
3633 end subroutine emis_read_climatology
3634
3635 !--------------------------------------------------------------------------
3636 ! emis_getIrEmissivity
3637 !--------------------------------------------------------------------------
3638 subroutine emis_getIrEmissivity (surfem1, nchn, sensorIndex, nprf, nchannels_max, sensorTovsIndexes)
3639 !
3640 !:Purpose: Assign new ir surface emissivities based on
3641 ! cmc analysis surface albedo, sea ice fraction and snow mask
3642 ! in addition to ceres surface type and water fraction.
3643 ! This is a subroutine that can apply to any instrument.
3644 !
3645 implicit none
3646
3647 ! Arguments:
3648 real(8), intent(out) :: surfem1(nchannels_max) ! IR surface emissivity estimate (0-1)
3649 integer, intent(in) :: nchn ! Number of channels
3650 integer, intent(in) :: sensorindex ! Sensor number
3651 integer, intent(in) :: nprf ! Number of profiles
3652 integer, intent(in) :: nchannels_max ! Total number of observations treated
3653 integer, intent(in) :: sensorTovsIndexes( nprf ) ! Profile position number
3654
3655 ! Locals:
3656 integer :: jc,jn
3657 integer :: ilat(nprf), ilon(nprf)
3658 real(8) :: latitudes(nprf), longitudes(nprf), satzang(nprf)
3659 real(8) :: wind_sfc(nprf), f_low(nprf), waven(nchn), em_oc(nchn,nprf), emi_mat(nchn,20)
3660
3661 ! Information to extract (transvidage)
3662 ! latitudes(nprf) -- latitude (-90 to 90)
3663 ! longitudes(nprf) -- longitude (0 to 360)
3664 ! satzang(nprf) -- satellite zenith angle (deg)
3665
3666 do jn = 1, nprf
3667 latitudes(jn) = tvs_profiles_nl(sensorTovsIndexes(jn)) % latitude
3668 longitudes(jn) = tvs_profiles_nl(sensorTovsIndexes(jn)) % longitude
3669 satzang(jn) = tvs_profiles_nl(sensorTovsIndexes(jn)) % zenangle
3670 end do
3671
3672 ! Assign surface properties from grid to profiles
3673 call interp_sfc(ilat,ilon, nprf,latitudes,longitudes,sensorTovsIndexes)
3674
3675
3676 ! Find the sensor bands (central) wavenumbers
3677 do jc = 1, nchn
3678 waven(jc) = tvs_coefs(sensorIndex) % coef % ff_cwn(jc)
3679 end do
3680
3681
3682 ! Get the CERES emissivity matrix for all sensor wavenumbers and surface types
3683 call ceres_ematrix(emi_mat, waven,nchn)
3684
3685
3686 ! Refine water emissivities
3687
3688 do jn = 1, nprf
3689 ! find surface wind
3690 wind_sfc(jn) = min(sqrt(tvs_profiles_nl(sensorTovsIndexes(jn)) % s2m %u**2 + tvs_profiles_nl(sensorTovsIndexes(jn)) % s2m % v**2 + 1.d-12),15.d0)
3691 end do
3692
3693 ! find new ocean emissivities
3694
3695 do jc = 1, nchn
3696 em_oc(jc,:)= emi_mat(jc,17)
3697 end do
3698
3699 call emi_sea (em_oc, waven,satzang,wind_sfc,nprf,nchn)
3700
3701
3702 ! Get surface emissivities
3703
3704 do jn = 1, nprf
3705 ! set albedo to 0.6 where snow is present
3706 if ( tvs_profiles_nl(sensorTovsIndexes(jn)) % skin % surftype == surftype_land .and. tvs_surfaceParameters(sensorTovsIndexes(jn)) % snow > 0.999 ) tvs_surfaceParameters(sensorTovsIndexes(jn)) % albedo = 0.6
3707 ! if albedo too high no water
3708 if ( tvs_surfaceParameters(sensorTovsIndexes(jn)) % albedo >= 0.55 ) tvs_surfaceParameters(sensorTovsIndexes(jn)) % pcnt_wat = 0.
3709 ! if water and CMC ice present then sea ice
3710 if ( tvs_profiles_nl(sensorTovsIndexes(jn)) % skin % surftype == surftype_sea .and. tvs_surfaceParameters(sensorTovsIndexes(jn)) % ice > 0.001 ) tvs_surfaceParameters(sensorTovsIndexes(jn)) % ltype = 20
3711 ! if land and CMC snow present then snow
3712 if ( tvs_profiles_nl(sensorTovsIndexes(jn)) % skin % surftype == surftype_land .and. tvs_surfaceParameters(sensorTovsIndexes(jn)) % snow > 0.999 ) tvs_surfaceParameters(sensorTovsIndexes(jn)) % ltype = 15
3713 do jc=1,nchn
3714 surfem1((jn-1)*nchn+jc) = tvs_surfaceParameters(sensorTovsIndexes(jn)) % pcnt_wat * em_oc(jc,jn) + &
3715 ( 1.d0 - tvs_surfaceParameters(sensorTovsIndexes(jn)) % pcnt_wat ) * emi_mat(jc,tvs_surfaceParameters(sensorTovsIndexes(jn)) % ltype)
3716 end do
3717 end do
3718
3719 ! Find the regional water fraction (here in a 15x15 pixel box centered on profile)
3720 call pcnt_box (f_low, waterFraction,nprf,ilat,ilon,kslat,kslon,7)
3721
3722 do jn = 1, nprf
3723 tvs_surfaceParameters(sensorTovsIndexes(jn)) % pcnt_reg = f_low(jn)
3724 end do
3725
3726 end subroutine emis_getIrEmissivity
3727
3728 !--------------------------------------------------------------------------
3729 ! interp_sfc
3730 !--------------------------------------------------------------------------
3731 subroutine interp_sfc (ilat, ilon, nprf, latitudes, longitudes, sensorTovsIndexes)
3732 !
3733 !:Purpose: Associate surface albedo, ice fraction, snow depth
3734 ! and ceres surface type and water fraction to observations profiles.
3735 !
3736 implicit none
3737
3738 ! Arguments:
3739 integer, intent(out) :: ilat(nprf) ! y-coordinate of profile
3740 integer, intent(out) :: ilon(nprf) ! x-coordinate of profile
3741 integer, intent(in) :: nprf ! number of profiles
3742 real(8), intent(in) :: latitudes(nprf) ! latitude (-90s to 90n)
3743 real(8), intent(in) :: longitudes(nprf) ! longitude (0 to 360)
3744 integer, intent(in) :: sensorTovsIndexes(nprf) ! observation index
3745
3746 ! Locals:
3747 character(len=20) :: cfile3,cfile5
3748 integer :: iun3,iun5
3749 integer :: iv7
3750 integer :: ix1,ix2,ix3,ix4,ix5, ix8,ix9,ix10,ix11,ix12
3751 integer :: iy3,iy4,iy5, iy8,iy9,iy10
3752 integer :: iz1,iz2,iz3,iz4,iz5, iz8,iz9,iz10,iz11,iz12
3753 integer :: ni3,nj3,nk3
3754 integer :: ni4,nj4,nk4
3755 integer :: ni5,nj5,nk5
3756 integer :: dateo,deet,npas,nbits,datyp
3757 integer :: ip1,ip2,ip3
3758 integer :: ig13,ig23,ig33,ig43
3759 integer :: ig14,ig24,ig34,ig44
3760 integer :: ig15,ig25,ig35,ig45
3761 integer :: swa,lng,dltf,ubc,ex1,ex2,ex3
3762 integer :: jn
3763 character(len=1) :: typvar
3764 character(len=1) :: grtyp3,grtyp4,grtyp5
3765 character(len=2) :: nomvar, snowvar
3766 character(len=8) :: etiket
3767 integer, external :: fnom,fstouv,fstinf,fstprm,fstfrm,fclos
3768 integer, external :: ezqkdef,ezdefset
3769 real(8) :: zig1,zig2,zig3,zig4
3770 integer :: ig1obs,ig2obs,ig3obs,ig4obs
3771 real (8) :: alat, alon, zzlat, zzlon
3772 ! fields on input grid
3773 real(8), allocatable :: glace(:,:), neige(:,:), alb(:,:)
3774 ! fields on output grid
3775 real(8) :: glace_intrpl(nprf,1), neige_intrpl(nprf,1), alb_intrpl(nprf,1)
3776
3777 ! printout header
3778 write(*,*)
3779 write(*,*) 'SUBROUTINE interp_sfc'
3780 write(*,*) '---------------------'
3781 write(*,*) ' called multiple time by bunch of ',nprf,' profiles'
3782 write(*,*) ' <RETURN CODES> SHOULD NOT BE NEGATIVE'
3783 write(*,*) '---------------------------------------------------'
3784
3785
3786 ! --- FOR CERES VARIABLES -------------
3787 ! Get number of pixels per degree of lat or lon
3788
3789 alat = dble(kslat)/180.d0
3790 alon = dble(kslon)/360.d0
3791
3792 do jn=1, nprf
3793
3794 ! get lat and lon within limits if necessary
3795 zzlat = min(latitudes(jn),89.999d0)
3796 zzlat = max(Zzlat,-89.999d0)
3797
3798 zzlon = min(longitudes(jn),359.999d0)
3799 zzlon = max(zzlon,0.d0)
3800
3801 ! Find in which surface field pixel is located the observation profile
3802
3803 ! Note : CERES grid at 1/6 resolution
3804 ! N-S : starts at N pole and excludes S pole
3805 ! W-E : starts at longitude 0 and excludes longitude 360
3806
3807 ilat(jn) = max( nint((zzlat + 90.d0) * alat),1)
3808 ilon(jn) = nint(zzlon * alon) + 1
3809 if (ilon(jn) > kslon) ilon(jn) = 1
3810
3811 ! Assign surface caracteristics to observation profiles
3812
3813 tvs_surfaceParameters(sensorTovsIndexes(jn)) % ltype = surfaceType(ilon(jn),ilat(jn))
3814 tvs_surfaceParameters(sensorTovsIndexes(jn)) % pcnt_wat = waterFraction(ilon(jn),ilat(jn))
3815
3816 end do
3817
3818 ! For ice, snow and albedo variables -------------
3819
3820 iun3 = 0
3821 iun5 = 0
3822
3823 ! File names
3824 cfile3 = 'sfc4airs' ! for ice fraction and snow cover
3825 cfile5 = 'sfc4airs_newalb' ! for albedo
3826
3827
3828 ! fnom: make the connections with the external files name
3829 ! success = 0
3830 write(*,*)
3831 ix1 = fnom(iun3,cfile3,'RND+R/O',0)
3832 write(*,*) 'file = sfc4airs : fnom : return = ', ix1
3833
3834 iz1 = fnom(iun5,cfile5,'RND+R/O',0)
3835 write(*,*) 'file = sfc4airs_newalb : fnom : return = ', iz1
3836
3837
3838 ! fstouv: open the standard files
3839 ! success = number of records found in the file
3840 write(*,*)
3841 ix2 = fstouv(iun3,'RND')
3842 write(*,*) 'file = sfc4airs : fstouv : return = ', ix2
3843 iz2 = fstouv(iun5,'RND')
3844 write(*,*) 'file = sfc4airs_newalb : fstouv : return = ', iz2
3845
3846
3847 ! fstinf: locate the records that matches the search keys
3848 ! success = handle of the record found after the search
3849 ! desired output = handle
3850 write(*,*)
3851 ix3 = fstinf(iun3,ni3,nj3,nk3,-1,'',-1,-1,-1,'','LG')
3852 write(*,*) 'variable = LG : fstinf : return = ', ix3
3853
3854 snowvar='SD'
3855 iy3 = fstinf(iun3,ni4,nj4,nk4,-1,'',-1,-1,-1,'',snowvar)
3856 write(*,*) 'variable = ', snowvar, ' : fstinf : return = ', iy3
3857 if ( iy3 < 0 ) then
3858 write(*,*) 'did not find ''SD'' so look for ''NE'''
3859 snowvar='NE'
3860 iy3 = fstinf(iun3,ni4,nj4,nk4,-1,'',-1,-1,-1,'',snowvar)
3861 write(*,*) 'variable = ', snowvar, ' : fstinf : return = ', iy3
3862 end if
3863
3864 iz3 = fstinf(iun5,ni5,nj5,nk5,-1,'',-1,-1,-1,'','AL')
3865 write(*,*) 'variable = AL : fstinf : return = ', iz3
3866
3867
3868 ! fstprm: get the description informations of the record given the key
3869 ! success = 0
3870 ! desired output = nix,njx,grtypx,igxx,ig1x,ig2x,ig3x,ig4x
3871
3872 write(*,*)
3873 ix4 = fstprm(ix3, dateo,deet,npas,ni3,nj3,nk3,nbits,datyp, &
3874 ip1,ip2,ip3,typvar,nomvar,etiket,grtyp3, &
3875 ig13,ig23,ig33,ig43,swa,lng,dltf,ubc,ex1,ex2,ex3)
3876 write(*,*) 'variable = LG : fstprm : return = ', ix4
3877
3878 iy4 = fstprm(iy3, dateo,deet,npas,ni4,nj4,nk4,nbits,datyp, &
3879 ip1,ip2,ip3,typvar,nomvar,etiket,grtyp4, &
3880 ig14,ig24,ig34,ig44,swa,lng,dltf,ubc,ex1,ex2,ex3)
3881 write(*,*) 'variable = ', snowvar, ' : fstprm : return = ', iy4
3882
3883 iz4 = fstprm(iz3, dateo,deet,npas,ni5,nj5,nk5,nbits,datyp, &
3884 ip1,ip2,ip3,typvar,nomvar,etiket,grtyp5, &
3885 ig15,ig25,ig35,ig45,swa,lng,dltf,ubc,ex1,ex2,ex3)
3886 write(*,*) 'variable = AL : fstprm : return = ', iz4
3887
3888
3889 ! allocation of the field on the grid
3890 allocate ( glace (ni3,nj3) )
3891 allocate ( neige (ni4,nj4) )
3892 allocate ( alb (ni5,nj5) )
3893
3894
3895 ! utl_fstlir: read records data (field on the grid) given the key
3896 ! success = handle of the record
3897 ! desired output = FIELD
3898 write(*,*)
3899
3900 ix5 = utl_fstlir(glace, iun3,ni3,nj3,nk3,-1,'',-1,-1,-1,'','LG')
3901 write(*,*) 'variable = LG : utl_fstlir : return = ', ix5
3902 iy5 = utl_fstlir(neige, iun3,ni4,nj4,nk4,-1,'',-1,-1,-1,'',snowvar)
3903 write(*,*) 'variable = ', snowvar, ' : utl_fstlir : return = ', iy5
3904 iz5 = utl_fstlir(alb, iun5,ni5,nj5,nk5,-1,'',-1,-1,-1,'','AL')
3905 write(*,*) 'variable = AL : utl_fstlir : return = ', iz5
3906
3907 ! int_CXGAIG: define the grid descriptors (integer form) of the
3908 ! observation profile output grid
3909 ! desired output = ig1OBS, ig2OBS, ig3OBS, ig4OBS
3910 zig1 = 0.0d0
3911 zig2 = 0.0d0
3912 zig3 = 1.0d0
3913 zig4 = 1.0d0
3914
3915 call int_cxgaig('L',ig1OBS,ig2OBS,ig3OBS,ig4OBS,zig1,zig2,zig3,zig4)
3916
3917
3918 ! int_EZGDEF: define the grid of the observations profiles (output grid)
3919 ! of type Y containing the lat-lon of profiles
3920 ! success = token to identify the grid
3921 ! desired output = token
3922 write(*,*)
3923 iv7 = int_ezgdef(nprf,1,'Y','L',ig1obs,ig2obs,ig3obs,ig4obs,longitudes,latitudes)
3924 write(*,*) 'apply to all variables : int_EZGDEF : return = ', iv7
3925
3926
3927 ! EZQKDEF: define the grid of the records data (input grid)
3928 ! success = token to identify the grid
3929 ! desired output = token
3930 ! EZDEFSET: interpolate from input grids to output grid
3931 ! success = key
3932 ! int_hInterpScalar: interpolation of the field on the input grid to observation profiles
3933 ! success = 0
3934 ! desired output = FIELD_intrpl
3935 write(*,*)
3936 ix8 = ezqkdef(ni3,nj3,grtyp3,ig13,ig23,ig33,ig43,iun3)
3937 write(*,*) 'variable = LG : ezqkdef : return = ', ix8
3938
3939 ix9 = ezdefset(iv7,ix8)
3940 write(*,*) 'variable = LG : ezdefset : return = ', ix9
3941
3942 ix10 = int_hInterpScalar(glace_intrpl,glace,interpDegree='NEAREST')
3943 write(*,*) 'variable = LG : int_hInterpScalar : return = ', ix10
3944
3945 write(*,*)
3946
3947 iy8 = ezqkdef(ni4,nj4,grtyp4,ig14,ig24,ig34,ig44,iun3)
3948 write(*,*) 'variable = ', snowvar, ' : ezqkdef : return = ', iy8
3949
3950 iy9 = ezdefset(iv7,iy8)
3951 write(*,*) 'variable = ', snowvar, ' : ezdefset : return = ', iy9
3952
3953 iy10 = int_hInterpScalar(neige_intrpl,neige,interpDegree='NEAREST')
3954 write(*,*) 'variable = ', snowvar, ' : int_hInterpScalar : return = ', iy10
3955
3956 write(*,*)
3957
3958 iz8 = ezqkdef(ni5,nj5,grtyp5,ig15,ig25,ig35,ig45,iun5)
3959 write(*,*) 'variable = AL : ezqkdef : return = ', iz8
3960
3961 iz9 = ezdefset(iv7,iz8)
3962 write(*,*) 'variable = AL : ezdefset : return = ', iz9
3963
3964 iz10 = int_hInterpScalar(alb_intrpl,alb,interpDegree='NEAREST')
3965 write(*,*) 'variable = AL : int_hInterpScalar : return = ', iz10
3966
3967
3968 ! fstfrm: close the standard files
3969 ! success = 0
3970 write(*,*)
3971 ix11 = fstfrm(iun3)
3972 write(*,*) 'file = sfc4airs : fstfrm : return = ', ix11
3973
3974 iz11 = fstfrm(iun5)
3975 write(*,*) 'file = sfc4airs_newalb : fstfrm : return = ', iz11
3976
3977
3978 ! fclos: release the connections with the external files name
3979 ! success = 0
3980
3981 write(*,*)
3982
3983 ix12 = fclos(iun3)
3984 write(*,*) 'file = sfc4airs : fclos : return = ', ix12
3985
3986 iz12 = fclos(iun5)
3987 write(*,*) 'file = sfc4airs_newalb : fclos : return = ', iz12
3988
3989 ! assign surface caracteristics to observation profiles
3990
3991 do jn=1, nprf
3992 tvs_surfaceParameters(sensorTovsIndexes(jn)) % ice = glace_intrpl(jn,1)
3993 tvs_surfaceParameters(sensorTovsIndexes(jn)) % snow = neige_intrpl(jn,1)
3994 tvs_surfaceParameters(sensorTovsIndexes(jn)) % albedo = alb_intrpl(jn,1)
3995 end do
3996
3997 deallocate(glace,neige,alb)
3998
3999 end subroutine interp_sfc
4000
4001 !--------------------------------------------------------------------------
4002 ! ceres_ematrix
4003 !--------------------------------------------------------------------------
4004 subroutine ceres_ematrix(emi_mat, waven, nchn)
4005 !
4006 !:Purpose: Set up emissivity versus fixed wavenumbers and surface types.
4007 !
4008 !:CERES:
4009 ! Emissivity data available at low spectral resolution: only 14 values
4010 ! to cover the entire spectrum. Thus, this can be used as a nominal value.
4011 ! The error associated with this emissivity can roughly be estimated to
4012 ! increase with lower emissivity as : (1-EMI)*0.5
4013 ! (i.e. as large as 0.10 for EMI=0.80 but better than 0.01 for EMI > 0.98)
4014 ! -No dependence on viewing angle is assumed.
4015 ! -Not to be used for oceans uncovered by ice.
4016 !
4017 !:Longwave Emmissivities in 12 original Fu bands + 2 extra to cover the range:
4018 !
4019 ! Longwave spectral intervals [cm-1] for the Fu & Liou code.
4020 !
4021 ! ==== ========== ========== ========== =========== ========== ========== ========= ========= ========= ========= ========= =============
4022 ! Band 1 2 3 4 5 6 7 8 9 10 11 12
4023 ! 2200-1900 1900-1700 1700-1400 1400-1250 1250-1100 1100-980 980-800 800-670 670-540 540-400 400-280 280-0
4024 ! ==== ========== ========== ========== =========== ========== ========== ========= ========= ========= ========= ========= =============
4025 !
4026 ! Two additional LW spectral intervals have been added in beyond 2200cm-1.
4027 !
4028 ! ===== =========== ===========
4029 ! Band 13 14
4030 ! 2500-2200 2850-2500
4031 ! ===== =========== ===========
4032 !
4033 ! Emissivity ems(band(1)) from April data, Table2 of Chen et al
4034 ! 11th Conf Sat Met, Madison, WI, p 514
4035 ! here regoganized as 14 13 1 2 ... 12 above
4036 !
4037 !:20 surface types:
4038 !
4039 ! =================== =================== =================== =====================
4040 ! 1= evergreen nleaf 2= evergreen bleaf 3= deciduous nleaf 4= deciduous bleaf
4041 ! 5= mixed forests 6= closed shrubs 7= open shrubs 8= woody savanna
4042 ! 9= savanna 10= grasslands 11= perma wet 12= croplands
4043 ! 13= urban 14= mosaic 15= snow 16= barren (deserts)
4044 ! 17= water 18= toundra 19= fresh snow 20= sea ice
4045 ! =================== =================== =================== =====================
4046 !
4047 implicit none
4048
4049 ! Arguments:
4050 integer, intent(in) :: nchn ! number of bands for which emissivity is needed
4051 real(8), intent(in) :: waven(nchn) ! wavenumbers (cm-1)
4052 real(8), intent(out) :: emi_mat(nchn, 20) ! emissivity (0.0-1.0)
4053
4054 ! Locals:
4055 integer :: i, nc, nt
4056 real(8) :: dum
4057
4058 ! CERES bands central wavenumber (covers 3.7 micron to 71.4 mic)
4059 integer, parameter :: nb=14
4060 real(8), parameter :: mid(nb) =(/ &
4061 2675.d0, 2350.d0, 2050.d0, 1800.d0, 1550.d0, 1325.d0, 1175.d0, 1040.d0, &
4062 890.d0, 735.d0, 605.d0, 470.d0, 340.d0, 140.d0 /)
4063
4064 ! CERES emissivity per wavenumber and surface types
4065 real(8), parameter :: emi_tab(nb,20)=(/ &
4066 0.951d0, 0.989d0, 0.989d0, 0.989d0, 0.990d0, 0.991d0, 0.991d0, 0.990d0, &
4067 0.990d0, 0.995d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4068 0.956d0, 0.989d0, 0.989d0, 0.989d0, 0.990d0, 0.991d0, 0.991d0, 0.990d0, &
4069 0.990d0, 0.995d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4070 0.929d0, 0.985d0, 0.985d0, 0.986d0, 0.984d0, 0.983d0, 0.979d0, 0.980d0, &
4071 0.973d0, 0.987d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4072 0.943d0, 0.985d0, 0.985d0, 0.986d0, 0.984d0, 0.983d0, 0.979d0, 0.980d0, &
4073 0.973d0, 0.987d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4074 0.945d0, 0.987d0, 0.987d0, 0.987d0, 0.987d0, 0.987d0, 0.985d0, 0.985d0, &
4075 0.982d0, 0.991d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4076 0.933d0, 0.949d0, 0.949d0, 0.970d0, 0.974d0, 0.971d0, 0.947d0, 0.958d0, &
4077 0.966d0, 0.975d0, 0.984d0, 0.984d0, 0.984d0, 0.984d0, &
4078 0.873d0, 0.873d0, 0.873d0, 0.934d0, 0.944d0, 0.939d0, 0.873d0, 0.904d0, &
4079 0.936d0, 0.942d0, 0.951d0, 0.951d0, 0.951d0, 0.951d0, &
4080 0.930d0, 0.987d0, 0.987d0, 0.990d0, 0.992d0, 0.993d0, 0.983d0, 0.975d0, &
4081 0.985d0, 0.993d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4082 0.926d0, 0.987d0, 0.987d0, 0.990d0, 0.992d0, 0.993d0, 0.983d0, 0.975d0, &
4083 0.985d0, 0.993d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4084 0.899d0, 0.987d0, 0.987d0, 0.990d0, 0.992d0, 0.993d0, 0.983d0, 0.975d0, &
4085 0.985d0, 0.993d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4086 0.951d0, 0.983d0, 0.983d0, 0.987d0, 0.987d0, 0.988d0, 0.983d0, 0.981d0, &
4087 0.987d0, 0.982d0, 0.986d0, 0.986d0, 0.986d0, 0.986d0, &
4088 0.924d0, 0.987d0, 0.987d0, 0.990d0, 0.992d0, 0.993d0, 0.983d0, 0.975d0, &
4089 0.985d0, 0.993d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4090 0.929d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4091 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4092 0.926d0, 0.987d0, 0.987d0, 0.989d0, 0.989d0, 0.990d0, 0.984d0, 0.980d0, &
4093 0.983d0, 0.992d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4094 0.972d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, 1.000d0, &
4095 1.000d0, 0.999d0, 0.999d0, 0.999d0, 0.999d0, 0.999d0, &
4096 0.866d0, 0.835d0, 0.835d0, 0.916d0, 0.934d0, 0.923d0, 0.835d0, 0.877d0, &
4097 0.921d0, 0.926d0, 0.934d0, 0.934d0, 0.934d0, 0.934d0, &
4098 0.973d0, 0.979d0, 0.979d0, 0.983d0, 0.982d0, 0.982d0, 0.984d0, 0.987d0, &
4099 0.989d0, 0.972d0, 0.972d0, 0.972d0, 0.972d0, 0.972d0, &
4100 0.968d0, 0.947d0, 0.947d0, 0.967d0, 0.988d0, 0.979d0, 0.975d0, 0.977d0, &
4101 0.992d0, 0.989d0, 0.989d0, 0.989d0, 0.989d0, 0.989d0, &
4102 0.984d0, 0.988d0, 0.988d0, 0.988d0, 0.988d0, 0.988d0, 0.988d0, 0.988d0, &
4103 0.988d0, 0.988d0, 0.988d0, 0.988d0, 0.988d0, 0.988d0, &
4104 0.964d0, 0.979d0, 0.979d0, 0.979d0, 0.979d0, 0.979d0, 0.979d0, 0.979d0, &
4105 0.979d0, 0.979d0, 0.979d0, 0.979d0, 0.979d0, 0.979d0 /)
4106
4107 do nt = 1, 20
4108 do nc = 1, nchn
4109 if ( waven(nc) > mid(1) ) then
4110 emi_mat(nc,nt) = emi_tab(1,nt)
4111 else if ( waven(nc) < mid(nb) ) then
4112 emi_mat(nc,nt) = emi_tab(nb,nt)
4113 else
4114 do i = 1, nb - 1
4115 if ( waven(nc) <= mid(i) .and. waven(nc) >= mid(i + 1) ) then
4116 dum = ( waven(nc) - mid(i) ) / ( mid(i + 1) - mid(i) )
4117 emi_mat(nc,nt) = emi_tab(i,nt) + ( emi_tab(i + 1,nt) - emi_tab(i,nt) ) * dum
4118 exit
4119 end if
4120 end do
4121 end if
4122 end do
4123 end do
4124
4125
4126 end subroutine ceres_ematrix
4127
4128 !--------------------------------------------------------------------------
4129 ! emi_sea
4130 !--------------------------------------------------------------------------
4131 subroutine emi_sea(em_oc, wnum, angle, wind, np, nc)
4132 !
4133 !:Purpose: GET OCEAN SURFACE EMISSIVITY
4134 !:Note: IMEM(NC), set to zero initially, on next call IMEM will have the
4135 ! right boundary channel to save search time in interpolation.
4136 ! IOPT=1 means activate IMEM option (all calls ask for same channels)
4137 !
4138 ! To get surface ocean emissivity for a group of channels with
4139 ! wavenumbers WNUM (cm-1) looking at one point with surface
4140 ! wind speed wind from angle angle.
4141 ! Based on Masuda,1988, Remote Sens. of Envir, 313-329.
4142 ! Coded emissivity routine based on Masuda's data by Tom Kleespies
4143 ! Covers 650-2857 cm-1 or 3.1-15.4 microns
4144 !
4145 !:CAUTION: extrapolated values from 769-650 cm-1
4146 ! and interpolated values between 2439-1250 cm-1
4147 !
4148 implicit none
4149
4150 ! Arguments:
4151 real(8), intent(out) :: em_oc(nc,np) ! Ocean emissivities (0.-1.)
4152 real(8), intent(in) :: wnum(nc) ! Channel wavenumbers (cm-1)
4153 real(8), intent(in) :: angle(np) ! Viewing angle (deg)
4154 real(8), intent(in) :: wind(np) ! Surface wind speed (m/s)
4155 integer, intent(in) :: np ! Number of profiles
4156 integer, intent(in) :: nc ! Number of channels
4157
4158 ! Locals:
4159 integer :: i, k, l
4160 integer :: imem(nc)
4161 integer :: mchan(2)
4162 real(8) :: dum
4163 real(8) :: emi2(2,np)
4164
4165 ! Masuda's 19 wavelengths converted to wavenumber
4166 real(8), parameter :: refw(19)=(/ 2857.1d0, 2777.7d0, 2702.7d0, 2631.6d0, 2564.1d0, &
4167 2500.0d0, 2439.0d0, 1250.0d0, 1190.5d0, 1136.3d0, &
4168 1087.0d0, 1041.7d0, 1000.0d0, 952.38d0, 909.09d0, &
4169 869.57d0, 833.33d0, 800.00d0, 769.23d0/)
4170
4171 ! imem options
4172 imem(:) = 0
4173
4174 do I = 1, nc
4175
4176 ! out of range
4177 if ( wnum(I) < 645.d0 .or. wnum(I) > refw(1) ) then
4178 write(*,'(A,1x,e12.4)') ' fatal: wavenumber out of range in emi_sea', wnum(I)
4179 stop
4180 else if ( wnum(I) <= refw(19) .and. wnum(I) > 645.d0 ) then
4181 ! extrapolated from 769 cm-1 to 645 cm-1: NOT FROM REAL DATA
4182 ! nevertheless thought to be much better than unity
4183 ! this is a region of relatively rapid emissivity change
4184 ! worst estimates for 700-645 cm-1, but these channels do not
4185 ! see the surface (strong co2 absorption).
4186 imem(I) = 18
4187 else
4188 ! CAUTION interpolation on large interval 1250-2439 cm-1
4189 ! where no data is available except that of ASTER. ASTER
4190 ! shows a relatively smooth variation with wavelength except
4191 ! for a sharp drop at 1600 cm-1 with highs at 1550 and 1650 cm-1
4192 ! with peak-to-peak variation of 1.5% in that narrow range.
4193 ! Worst estimates would be between 1400-1800 cm-1 in HIRS ch 12
4194 ! which only in very cold atmospheres sees the surface.
4195 do k = 1, 18
4196 if ( wnum(I) > refw(k + 1) .and. wnum(I) <= refw(k) ) then
4197 imem(I) = k
4198 end if
4199 end do
4200
4201 end if
4202
4203 mchan(1)= imem(I)
4204 mchan(2)= imem(I) + 1
4205
4206 dum = ( wnum(I) - refw(mchan(1)) ) / ( refw(mchan(2)) - refw(mchan(1)) )
4207
4208 call COMP_IR_EMISS(emi2, wind,angle,2,np,mchan)
4209
4210 ! interpolation/extrapolation in wavenumber
4211
4212 do L = 1, np
4213
4214 em_oc(I,L) = emi2(1,L) + ( emi2(2,L) - emi2(1,L) ) * dum
4215
4216 end do
4217
4218 end do
4219
4220
4221 end subroutine emi_sea
4222
4223
4224 !--------------------------------------------------------------------------
4225 ! tvs_getCommonChannelSet
4226 !--------------------------------------------------------------------------
4227 subroutine tvs_getCommonChannelSet(channels,countUniqueChannel, listAll)
4228 !
4229 !:Purpose: get common channels among all MPI tasks
4230 !
4231 implicit none
4232
4233 ! Arguments:
4234 integer, intent(in) :: channels(:)
4235 integer, intent(out) :: countUniqueChannel
4236 integer, intent(out) :: listAll(:)
4237
4238 ! Locals:
4239 integer :: channelsb(tvs_maxChannelNumber)
4240 integer :: ierr, i, j
4241 integer, allocatable :: listGlobal(:)
4242 logical :: found
4243
4244 if (size(channels) > tvs_maxChannelNumber) then
4245 write(*,*) 'You need to increase tvs_maxChannelNumber in tovsNL_mod !',size(channels), tvs_maxChannelNumber
4246 call utl_abort('tvs_getCommonChannelSet')
4247 end if
4248
4249 if (mmpi_myid ==0) then
4250 allocate(listGlobal(mmpi_nprocs*tvs_maxChannelNumber))
4251 else
4252 allocate(listGlobal(1))
4253 end if
4254
4255 listAll(:) = 0
4256 listGlobal(:) = 0
4257 channelsb(:) = 0
4258 channelsb(1:size(channels)) = channels(:)
4259
4260 call rpn_comm_barrier('GRID',ierr)
4261
4262 call rpn_comm_gather(channelsb, tvs_maxChannelNumber, 'MPI_INTEGER', listGlobal, &
4263 tvs_maxChannelNumber, 'MPI_INTEGER', 0, 'GRID', ierr)
4264 countUniqueChannel = 0
4265 if ( mmpi_myid == 0 ) then
4266 call isort(listGlobal, mmpi_nprocs*tvs_maxChannelNumber)
4267 do i=1, mmpi_nprocs * tvs_maxChannelNumber
4268 if (listGlobal(i) > 0) then
4269 found = .false.
4270 LOOPJ: do j=countUniqueChannel,1,-1
4271 if (listGlobal(i) == listAll(j) ) then
4272 found =.true.
4273 exit LOOPJ
4274 end if
4275 end do LOOPJ
4276 if (.not.found) then
4277 countUniqueChannel = countUniqueChannel + 1
4278 listAll(countUniqueChannel) = listGlobal(i)
4279 end if
4280 end if
4281 end do
4282 end if
4283
4284 call rpn_comm_bcast(countUniqueChannel, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4285 call rpn_comm_bcast(listAll(1:countUniqueChannel), countUniqueChannel, 'MPI_INTEGER', 0, 'GRID', ierr)
4286
4287 deallocate(listGlobal)
4288
4289 end subroutine tvs_getCommonChannelSet
4290
4291
4292 !--------------------------------------------------------------------------
4293 ! tvs_rttov_read_coefs
4294 !--------------------------------------------------------------------------
4295 subroutine tvs_rttov_read_coefs(errorStatus, coefs, opts, channels, instrument)
4296 !
4297 !:Purpose: MPI wrapper for rttov_read_coefs
4298 ! the coefficient files are read by MPI task 0
4299 ! and then broadcasted to the other tasks according to the selected
4300 ! channels. Argument channels is mandatory (it is optional in rttov_setup)
4301 ! optional argument channels_rec was removed (it is useful only in principal component mode)
4302 ! other optionnal arguments were removed :
4303 ! * form_coef, to specify format
4304 ! * form_scaer,
4305 ! * form_sccld,
4306 ! * form_pccoef,
4307 ! * file_coef, to specify filename
4308 ! * file_scaer,
4309 ! * file_sccld,
4310 ! * file_pccoef,
4311 ! * file_id_coef, to specify fortran unit number
4312 ! * file_id_scaer,
4313 ! * file_id_sccld,
4314 ! * file_id_pccoef,
4315 ! * path, to specify the path to look for coefficient files
4316 !
4317 ! if necessary these arguments could be added (ask S. Heilliette)
4318 ! also this subroutine will work only for clear sky radiance computations
4319 ! if somebody wants to do realistic cloud or aerosol affected radiance simulations
4320 ! some changes are needed. Ask me in that case. (S. Heilliette)
4321 ! It is implicitely assumed that the options are the same for all MPI tasks for a given instrument
4322 ! No check will be done (options for task 0 will be used for all tasks).
4323 ! Only differences in channel lists are accounted for.
4324 !
4325 implicit none
4326
4327 ! Arguments:
4328 integer(kind=jpim), intent(out) :: errorStatus ! Error status
4329 type(rttov_coefs), intent(out) :: coefs ! Rttov coefficient structure
4330 type(rttov_options), intent(in) :: opts ! Rttov option structure
4331 integer(kind=jpim), intent(in) :: channels(:) ! Channel list
4332 integer(kind=jpim), intent(in) :: instrument(3) ! Instrument vector
4333
4334 ! Locals:
4335 real(8), allocatable :: bigArray(:,:,:,:)
4336 integer :: i, j, ichan, ierr, countUniqueChannel, indexchan(size(channels)), listAll(tvs_maxChannelNumber)
4337 logical :: associated0
4338 integer :: nlte_count, nlte_start,isol,isat,nlte_file_nchan
4339 integer, allocatable :: nlte_chans(:)
4340
4341 write(*,*) 'tvs_rttov_read_coefs: Starting'
4342
4343 ! First step: we should determine a common set of channels among MPI tasks
4344 call tvs_getCommonChannelSet(channels,countUniqueChannel, listAll)
4345
4346 ! Second step: mpi task 0 will do the job
4347 if ( mmpi_myid == 0 ) then
4348 call rttov_read_coefs ( &
4349 errorStatus, &! out
4350 coefs, &
4351 opts, &
4352 instrument=instrument, &! in
4353 channels=listAll(1:countUniqueChannel) ) ! in option
4354 if (errorStatus /= errorStatus_success) then
4355 write(*,*) 'tvs_rttov_read_coefs: failure in rttov_read_coefs while reading RTTOV coefficient file', instrument, errorStatus
4356 call utl_abort('tvs_rttov_read_coefs')
4357 end if
4358 else
4359 call rttov_nullify_coef(coefs%coef)
4360 end if
4361
4362 ! Third step: common (i.e. independent from the channel list) parameters are simply broadcasted to other processors
4363 ! Scalar and fixed size arrays and strings first
4364 coefs%initialised = .true. ! Logical flag for initialization
4365 call rpn_comm_bcast(coefs%coef%id_platform, 1, 'MPI_INTEGER', 0, 'GRID', ierr)! OK
4366 call rpn_comm_bcast(coefs%coef%id_sat, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! OK
4367 call rpn_comm_bcast(coefs%coef%id_inst, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! OK
4368 call rpn_comm_bcast(coefs%coef%id_sensor, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! OK
4369 call rpn_comm_bcast(coefs%coef%id_comp_lvl, 1, 'MPI_INTEGER', 0, 'GRID', ierr)! OK
4370 call rpn_comm_bcast(coefs%coef%id_comp_pc, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! OK
4371
4372 call rpn_comm_bcast(coefs%coef%fmv_model_ver, 1, 'MPI_INTEGER', 0, 'GRID', ierr) !ok
4373 call rpn_comm_bcast(coefs%coef%fmv_chn, 1, 'MPI_INTEGER', 0, 'GRID', ierr) !ok
4374 call rpn_comm_bcast(coefs%coef%fmv_gas, 1, 'MPI_INTEGER', 0, 'GRID', ierr) !ok
4375 call rpn_comm_bcast(coefs%coef%fmv_ori_nchn, 1, 'MPI_INTEGER', 0, 'GRID', ierr) !ok
4376
4377 call rpn_comm_bcast(coefs%coef%nmixed, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for Mixed Gases
4378 call rpn_comm_bcast(coefs%coef%nwater, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for Water Vapour
4379 call rpn_comm_bcast(coefs%coef%nozone, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for Ozone
4380 call rpn_comm_bcast(coefs%coef%nwvcont, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for WV continuum
4381 call rpn_comm_bcast(coefs%coef%nco2, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for CO2
4382 call rpn_comm_bcast(coefs%coef%nn2o , 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for N2O
4383 call rpn_comm_bcast(coefs%coef%nco, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for CO
4384 call rpn_comm_bcast(coefs%coef%nch4, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for CH4
4385 call rpn_comm_bcast(coefs%coef%nso2, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of variables/predictors for SO2
4386
4387 call rpn_comm_bcast(coefs%coef%nlevels, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of levels(pres/absorber) same for all gases
4388 call rpn_comm_bcast(coefs%coef%nlayers, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of layers(pres/absorber) nlevels-1
4389 call rpn_comm_bcast(coefs%coef%pmc_nlay, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4390 call rpn_comm_bcast(coefs%coef%pmc_nvar, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4391 call rpn_comm_bcast(coefs%coef%IncZeeman, 1, 'MPI_LOGICAL', 0, 'GRID', ierr) ! Flag to include Zeeman effect for this sensor
4392 call rpn_comm_bcast(coefs%coef%solarcoef, 1, 'MPI_LOGICAL', 0, 'GRID', ierr) ! Flag to include solar reflection
4393 call rpn_comm_bcast(coefs%coef%nltecoef, 1, 'MPI_LOGICAL', 0, 'GRID', ierr) ! Flag to include nlte corrections
4394 call rpn_comm_bcast(coefs%coef%pmc_shift, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
4395
4396 call rpn_comm_bcast(coefs%coef%ncmixed, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for Mixed Gases
4397 call rpn_comm_bcast(coefs%coef%ncwater, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for Water Vapour
4398 call rpn_comm_bcast(coefs%coef%ncozone, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for Ozone
4399 call rpn_comm_bcast(coefs%coef%ncwvcont, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for WV continuum
4400 call rpn_comm_bcast(coefs%coef%ncco2, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for CO2
4401 call rpn_comm_bcast(coefs%coef%ncn2o, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for N2O
4402 call rpn_comm_bcast(coefs%coef%ncco , 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for CO
4403 call rpn_comm_bcast(coefs%coef%ncch4, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for CH4
4404 call rpn_comm_bcast(coefs%coef%ncso2, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for SO2
4405
4406 call rpn_comm_bcast(coefs%coef%nccmixed, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for Mixed Gases
4407 call rpn_comm_bcast(coefs%coef%nccwater, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for Water Vapour
4408 call rpn_comm_bcast(coefs%coef%nccozone, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for Ozone
4409 call rpn_comm_bcast(coefs%coef%nccwvcont, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for WV continuum
4410 call rpn_comm_bcast(coefs%coef%nccco2, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for CO2
4411 call rpn_comm_bcast(coefs%coef%nccn2o, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for N2O
4412 call rpn_comm_bcast(coefs%coef%nccco , 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for CO
4413 call rpn_comm_bcast(coefs%coef%nccch4, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for CH4
4414 call rpn_comm_bcast(coefs%coef%nccso2, 1, 'MPI_INTEGER', 0, 'GRID', ierr) ! number of coefficients for SO2
4415
4416 call rpn_comm_bcast(coefs%coef%ws_nomega, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4417
4418 call rpn_comm_bcast(coefs%coef%id_creation_date, 3, 'MPI_INTEGER', 0, 'GRID', ierr) ! OK
4419 call rpn_comm_bcastc(coefs%coef%id_creation, 80, 'MPI_CHARACTER', 0, 'GRID', ierr) ! OK
4420 call rpn_comm_bcastc(coefs%coef%id_Common_name, 32, 'MPI_CHARACTER', 0, 'GRID', ierr) !OK
4421 do i=1, 100
4422 call rpn_comm_bcastc(coefs%coef%line_by_line(i), 132, 'MPI_CHARACTER', 0, 'GRID', ierr) !ok
4423 call rpn_comm_bcastc(coefs%coef%readme_srf(i), 132, 'MPI_CHARACTER', 0, 'GRID', ierr) !ok
4424 end do
4425 call rpn_comm_bcastc(coefs%coef%fmv_model_def, 32, 'MPI_CHARACTER', 0, 'GRID', ierr) !OK
4426 call rpn_comm_bcast(coefs%coef%fc_planck_c1, 1, 'MPI_REAL8', 0, 'GRID', ierr) ! first radiation constant (mW/(m2*sr*cm-4)) !ok
4427 call rpn_comm_bcast(coefs%coef%fc_planck_c2, 1, 'MPI_REAL8', 0, 'GRID', ierr) !second radiation constant (mW/(m2*sr*cm-4)) !ok
4428 call rpn_comm_bcast(coefs%coef%fc_sat_height, 1, 'MPI_REAL8', 0, 'GRID', ierr)! satellite nominal altitude (km) !ok
4429
4430 call rpn_comm_bcast(coefs%coef%pmc_lengthcell, 1, 'MPI_REAL8', 0, 'GRID', ierr)
4431 call rpn_comm_bcast(coefs%coef%pmc_tempcell, 1, 'MPI_REAL8', 0, 'GRID', ierr)
4432 call rpn_comm_bcast(coefs%coef%pmc_betaplus1, 1, 'MPI_REAL8', 0, 'GRID', ierr)
4433 ! FASTEM section
4434 call rpn_comm_bcast(coefs%coef%ssirem_ver, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4435 call rpn_comm_bcast(coefs%coef%iremis_version, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4436 call rpn_comm_bcast(coefs%coef%iremis_ncoef, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4437 call rpn_comm_bcast(coefs%coef%iremis_angle0, 1, 'MPI_REAL8', 0, 'GRID', ierr)
4438 call rpn_comm_bcast(coefs%coef%iremis_tskin0, 1, 'MPI_REAL8', 0, 'GRID', ierr)
4439 call rpn_comm_bcast(coefs%coef%ratoe, 1, 'MPI_REAL8', 0, 'GRID', ierr)
4440 ! then variable size vectors
4441 ! this one must be done first because it is used to dimension other ones ....
4442 if (mmpi_myid > 0) allocate( coefs%coef%fmv_lvl(coefs%coef%fmv_gas))
4443 call rpn_comm_bcast(coefs%coef%fmv_lvl, coefs%coef%fmv_gas, 'MPI_INTEGER', 0, 'GRID', ierr)
4444 if (mmpi_myid > 0) then
4445 if (coefs%coef%nltecoef) allocate( coefs%coef%nlte_coef)
4446 allocate( coefs%coef%fmv_gas_id(coefs%coef%fmv_gas))
4447 allocate( coefs%coef%fmv_gas_pos(ngases_max)) !size is from rttov consts
4448 allocate( coefs%coef%fmv_var(coefs%coef%fmv_gas))
4449 allocate( coefs%coef%fmv_coe(coefs%coef%fmv_gas)) ! number of coefficients by gas (fmv_gas)
4450 allocate( coefs%coef%fmv_ncorr(coefs%coef%fmv_gas)) ! number of coefs by gas for correction term (fmv_gas) (v13 only)
4451 allocate( coefs%coef%ws_npoint(coefs%coef%ws_nomega) )
4452 allocate( coefs%coef%ws_k_omega(coefs%coef%ws_nomega) )
4453 allocate( coefs%coef%ref_prfl_p(coefs%coef%nlevels ) )
4454 allocate( coefs%coef%ref_prfl_t(coefs%coef%nlevels, coefs%coef%fmv_gas ) )
4455 allocate( coefs%coef%ref_prfl_mr(coefs%coef%nlevels, coefs%coef%fmv_gas ) )
4456 allocate( coefs%coef%bkg_prfl_mr(coefs%coef%nlevels, coefs%coef%fmv_gas ) )
4457 allocate( coefs%coef%lim_prfl_p(coefs%coef%nlevels ) )
4458 allocate( coefs%coef%lim_prfl_tmax(coefs%coef%fmv_lvl(gas_id_mixed) ) )
4459 allocate( coefs%coef%lim_prfl_tmin(coefs%coef%fmv_lvl(gas_id_mixed) ) )
4460 allocate( coefs%coef%lim_prfl_gmax(coefs%coef%fmv_lvl(gas_id_mixed), coefs%coef%fmv_gas ) )
4461 allocate( coefs%coef%lim_prfl_gmin(coefs%coef%fmv_lvl(gas_id_mixed), coefs%coef%fmv_gas ) )
4462 allocate( coefs%coef%env_prfl_tmax(coefs%coef%nlevels ) )
4463 allocate( coefs%coef%env_prfl_tmin(coefs%coef%nlevels ) )
4464 allocate( coefs%coef%env_prfl_gmax(coefs%coef%nlevels, coefs%coef%fmv_gas ) )
4465 allocate( coefs%coef%env_prfl_gmin(coefs%coef%nlevels, coefs%coef%fmv_gas ) )
4466 allocate( coefs%coef%dpp(0:coefs%coef%nlayers) )
4467 allocate( coefs%coef%dp(coefs%coef%nlayers) )
4468 allocate( coefs%coef%tstar(coefs%coef%nlayers) )
4469 allocate( coefs%coef%tstar_r(coefs%coef%nlayers) )
4470 allocate( coefs%coef%tstar_wsum_r(0:coefs%coef%nlayers) )
4471 if (coefs%coef%fmv_model_ver == 8) &
4472 allocate( coefs%coef%tstarmod_wsum_r(coefs%coef%nlayers) )
4473 if (coefs%coef%fmv_model_ver <= 9) &
4474 allocate( coefs%coef%tstar_uwsum_r(0:coefs%coef%nlayers) )
4475 allocate( coefs%coef%wstar(coefs%coef%nlayers) )
4476 allocate( coefs%coef%wstar_r(coefs%coef%nlayers) )
4477 allocate( coefs%coef%wstar_wsum_r(0:coefs%coef%nlayers) )
4478 allocate( coefs%coef%wtstar_wsum_r(0:coefs%coef%nlayers) )
4479 end if
4480
4481 call broadcastI41dArray( coefs%coef%fmv_gas_id )
4482 call broadcastI41dArray( coefs%coef%fmv_gas_pos )
4483 call broadcastI41dArray( coefs%coef%fmv_var )
4484 call broadcastI41dArray( coefs%coef%fmv_coe )
4485 call broadcastI41dArray( coefs%coef%fmv_ncorr )
4486 call broadcastR81dArray( coefs%coef%ws_npoint )
4487 call broadcastR81dArray( coefs%coef%ws_k_omega )
4488 call broadcastR81dArray( coefs%coef%ref_prfl_p )
4489 call broadcastR82dArray( coefs%coef%ref_prfl_t )
4490 call broadcastR82dArray( coefs%coef%ref_prfl_mr )
4491 call broadcastR82dArray( coefs%coef%bkg_prfl_mr )
4492 call broadcastR81dArray( coefs%coef%lim_prfl_p )
4493 call broadcastR81dArray( coefs%coef%lim_prfl_tmax )
4494 call broadcastR81dArray( coefs%coef%lim_prfl_tmin )
4495 call broadcastR82dArray( coefs%coef%lim_prfl_gmax )
4496 call broadcastR82dArray( coefs%coef%lim_prfl_gmin )
4497 call broadcastR81dArray( coefs%coef%env_prfl_tmax )
4498 call broadcastR81dArray( coefs%coef%env_prfl_tmin )
4499 call broadcastR82dArray( coefs%coef%env_prfl_gmax )
4500 call broadcastR82dArray( coefs%coef%env_prfl_gmin )
4501 call broadcastR81dArray( coefs%coef%dp )
4502 call broadcastR81dArray( coefs%coef%dpp )
4503 call broadcastR81dArray( coefs%coef%tstar )
4504 call broadcastR81dArray( coefs%coef%tstar_r )
4505 call broadcastR81dArray( coefs%coef%tstar_wsum_r )
4506 if (coefs%coef%fmv_model_ver == 8) &
4507 call broadcastR81dArray( coefs%coef%tstarmod_wsum_r )
4508 if (coefs%coef%fmv_model_ver <= 9) &
4509 call broadcastR81dArray( coefs%coef%tstar_uwsum_r )
4510 call broadcastR81dArray( coefs%coef%wstar )
4511 call broadcastR81dArray( coefs%coef%wstar_r )
4512 call broadcastR81dArray( coefs%coef%wstar_wsum_r )
4513 call broadcastR81dArray( coefs%coef%wtstar_wsum_r )
4514 if (coefs%coef%nozone > 0) then
4515 if (mmpi_myid > 0) then
4516 allocate( coefs%coef%to3star(coefs%coef%nlayers) )
4517 allocate( coefs%coef%ostar(coefs%coef%nlayers) )
4518 allocate( coefs%coef%to3star_r(coefs%coef%nlayers) )
4519 allocate( coefs%coef%ostar_r(coefs%coef%nlayers) )
4520 allocate( coefs%coef%ostar_wsum_r(0:coefs%coef%nlayers) )
4521 end if
4522 call rpn_comm_bcast(coefs%coef%to3star, size(coefs%coef%to3star) , 'MPI_REAL8', 0, 'GRID', ierr)
4523 call rpn_comm_bcast(coefs%coef%ostar, size(coefs%coef%ostar) , 'MPI_REAL8', 0, 'GRID', ierr)
4524 call rpn_comm_bcast(coefs%coef%to3star_r, size(coefs%coef%to3star_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4525 call rpn_comm_bcast(coefs%coef%ostar_r, size(coefs%coef%ostar_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4526 call rpn_comm_bcast(coefs%coef%ostar_wsum_r, size(coefs%coef%ostar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4527 end if
4528 if ( coefs%coef%nco2 > 0) then
4529 if (mmpi_myid>0) then
4530 allocate( coefs%coef%co2star(coefs%coef%nlayers) )
4531 allocate( coefs%coef%co2star_r(coefs%coef%nlayers) )
4532 allocate( coefs%coef%co2star_wsum_r(0:coefs%coef%nlayers) )
4533 end if
4534 call rpn_comm_bcast(coefs%coef%co2star, size(coefs%coef%co2star) , 'MPI_REAL8', 0, 'GRID', ierr)
4535 call rpn_comm_bcast(coefs%coef%co2star_r, size(coefs%coef%co2star_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4536 call rpn_comm_bcast(coefs%coef%co2star_wsum_r, size(coefs%coef%co2star_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4537 end if
4538 if ( coefs%coef%nn2o > 0) then
4539 if (mmpi_myid>0) then
4540 allocate( coefs%coef%n2ostar(coefs%coef%nlayers) )
4541 allocate( coefs%coef%n2ostar_r(coefs%coef%nlayers) )
4542 allocate( coefs%coef%n2ostar_wsum_r(0:coefs%coef%nlayers) )
4543 allocate( coefs%coef%n2otstar_wsum_r(0:coefs%coef%nlayers) )
4544 end if
4545 call rpn_comm_bcast(coefs%coef%n2ostar, size(coefs%coef%n2ostar) , 'MPI_REAL8', 0, 'GRID', ierr)
4546 call rpn_comm_bcast(coefs%coef%n2ostar_r, size(coefs%coef%n2ostar_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4547 call rpn_comm_bcast(coefs%coef%n2ostar_wsum_r, size(coefs%coef%n2ostar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4548 call rpn_comm_bcast(coefs%coef%n2otstar_wsum_r, size(coefs%coef%n2otstar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4549 end if
4550 if ( coefs%coef%nco > 0) then
4551 if (mmpi_myid>0) then
4552 allocate( coefs%coef%costar(coefs%coef%nlayers) )
4553 allocate( coefs%coef%costar_r(coefs%coef%nlayers) )
4554 allocate( coefs%coef%costar_wsum_r(0:coefs%coef%nlayers) )
4555 allocate( coefs%coef%cotstar_wsum_r(0:coefs%coef%nlayers) )
4556 end if
4557 call rpn_comm_bcast(coefs%coef%costar, size(coefs%coef%costar) , 'MPI_REAL8', 0, 'GRID', ierr)
4558 call rpn_comm_bcast(coefs%coef%costar_r, size(coefs%coef%costar_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4559 call rpn_comm_bcast(coefs%coef%costar_wsum_r, size(coefs%coef%costar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4560 call rpn_comm_bcast(coefs%coef%cotstar_wsum_r, size(coefs%coef%cotstar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4561 end if
4562 if ( coefs%coef%nch4 > 0) then
4563 if (mmpi_myid>0) then
4564 allocate( coefs%coef%ch4star(coefs%coef%nlayers) )
4565 allocate( coefs%coef%ch4star_r(coefs%coef%nlayers) )
4566 allocate( coefs%coef%ch4star_wsum_r(0:coefs%coef%nlayers) )
4567 allocate( coefs%coef%ch4tstar_wsum_r(coefs%coef%nlayers) )
4568 end if
4569 call rpn_comm_bcast(coefs%coef%ch4star, size(coefs%coef%ch4star) , 'MPI_REAL8', 0, 'GRID', ierr)
4570 call rpn_comm_bcast(coefs%coef%ch4star_r, size(coefs%coef%ch4star_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4571 call rpn_comm_bcast(coefs%coef%ch4star_wsum_r, size(coefs%coef%ch4star_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4572 call rpn_comm_bcast(coefs%coef%ch4tstar_wsum_r, size(coefs%coef%ch4tstar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4573 end if
4574 if (coefs%coef%nso2 > 0) then
4575 if (mmpi_myid>0) then
4576 allocate( coefs%coef%so2star(coefs%coef%nlayers) )
4577 allocate( coefs%coef%so2star_r(coefs%coef%nlayers) )
4578 allocate( coefs%coef%so2star_wsum_r(0:coefs%coef%nlayers) )
4579 allocate( coefs%coef%so2tstar_wsum_r(coefs%coef%nlayers) )
4580 end if
4581 call rpn_comm_bcast(coefs%coef%so2star, size(coefs%coef%so2star) , 'MPI_REAL8', 0, 'GRID', ierr)
4582 call rpn_comm_bcast(coefs%coef%so2star_r, size(coefs%coef%so2star_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4583 call rpn_comm_bcast(coefs%coef%so2star_wsum_r, size(coefs%coef%so2star_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4584 call rpn_comm_bcast(coefs%coef%so2tstar_wsum_r, size(coefs%coef%so2tstar_wsum_r) , 'MPI_REAL8', 0, 'GRID', ierr)
4585 end if
4586 ! Fourth step: channel dependent parameters are extracted according to the channel list and sent to each MPI task
4587 coefs%coef%fmv_chn = size( channels )
4588
4589 if (mmpi_myid == 0) deallocate ( coefs%coef%ff_ori_chn)
4590 allocate ( coefs%coef%ff_ori_chn(coefs%coef%fmv_chn) )
4591 coefs%coef%ff_ori_chn = channels
4592
4593 do i=1, coefs%coef%fmv_chn
4594 loopj2:do j=1, countUniqueChannel
4595 if (listAll(j) == channels(i)) then
4596 indexchan(i) = j
4597 exit loopj2
4598 end if
4599 end do loopj2
4600 end do
4601
4602 ! 1D arrays first
4603 call extractI41dArray(coefs%coef%pw_val_chn, countUniqueChannel,indexchan)
4604 call extractI41dArray(coefs%coef%ff_val_chn, countUniqueChannel,indexchan)
4605 call extractR81dArray(coefs%coef%ff_cwn, countUniqueChannel,indexchan)
4606 call extractR81dArray(coefs%coef%ff_bco, countUniqueChannel,indexchan)
4607 call extractR81dArray(coefs%coef%ff_bcs, countUniqueChannel,indexchan)
4608 call extractR81dArray(coefs%coef%ff_gam, countUniqueChannel,indexchan)
4609 call extractI41dArray(coefs%coef%tt_val_chn, countUniqueChannel,indexchan)
4610 call extractR81dArray(coefs%coef%tt_a0, countUniqueChannel,indexchan)
4611 call extractR81dArray(coefs%coef%tt_a1, countUniqueChannel,indexchan)
4612 call extractI41dArray(coefs%coef%ss_val_chn, countUniqueChannel,indexchan)
4613 call extractR81dArray(coefs%coef%ss_solar_spectrum, countUniqueChannel,indexchan)
4614 if ( coefs%coef%fmv_model_ver > 9) then
4615 call extractR81dArray(coefs%coef%ss_rayleigh_ext, countUniqueChannel,indexchan)
4616 end if
4617 call extractCmplx81dArray(coefs%coef%woc_waopc_ow, countUniqueChannel, indexchan)
4618 call extractCmplx81dArray(coefs%coef%woc_waopc_fw, countUniqueChannel, indexchan)
4619
4620 call extractI41dArray(coefs%coef%fastem_polar, countUniqueChannel,indexchan)
4621 call extractR81dArray(coefs%coef%pol_phi, countUniqueChannel,indexchan)
4622 call extractR81dArray(coefs%coef%ssirem_a0, countUniqueChannel,indexchan)
4623 call extractR81dArray(coefs%coef%ssirem_a1, countUniqueChannel,indexchan)
4624 call extractR81dArray(coefs%coef%ssirem_a2, countUniqueChannel,indexchan)
4625 call extractR81dArray(coefs%coef%ssirem_xzn1, countUniqueChannel,indexchan)
4626 call extractR81dArray(coefs%coef%ssirem_xzn2, countUniqueChannel,indexchan)
4627
4628 call extractR81dArray(coefs%coef%planck1, countUniqueChannel,indexchan)
4629 call extractR81dArray(coefs%coef%planck2, countUniqueChannel,indexchan)
4630 if (coefs%coef%id_sensor == sensor_id_mw .or. coefs%coef%id_sensor == sensor_id_po) &
4631 call extractR81dArray(coefs%coef%frequency_ghz, countUniqueChannel,indexchan)
4632
4633 call extractR81dArray(coefs%coef%pmc_pnominal, countUniqueChannel,indexchan)
4634 call extractR81dArray(coefs%coef%pmc_ppmc, countUniqueChannel,indexchan)
4635 call extractR81dArray(coefs%coef%pol_fac_v, countUniqueChannel,indexchan)
4636 call extractR81dArray(coefs%coef%pol_fac_h, countUniqueChannel,indexchan)
4637
4638 ! 2D arrays
4639 call extractR82dArray(coefs%coef%iremis_coef,coefs%coef%iremis_ncoef,countUniqueChannel,indexchan)
4640 ! 3D arrays
4641 call extractR83dArray(coefs%coef%pmc_coef,coefs%coef%pmc_nlay,countUniqueChannel,coefs%coef%pmc_nvar,indexchan)
4642
4643 ! then coefficients. It is more complicated with RTTOV12
4644 call dispatch_fast_coef(errorStatus, coefs%coef%thermal, coefs%coef%fmv_gas_id, coefs%coef%fmv_coe, coefs%coef%fmv_model_ver, &
4645 coefs%coef%nlayers, coefs%coef%fmv_gas)
4646 if (coefs%coef%fmv_model_ver > 9) THEN
4647 call dispatch_fast_coef(errorStatus, coefs%coef%thermal_corr, coefs%coef%fmv_gas_id, coefs%coef%fmv_ncorr, coefs%coef%fmv_model_ver, &
4648 coefs%coef%nlayers, coefs%coef%fmv_gas)
4649 end if
4650 if (coefs%coef%solarcoef) then
4651 call dispatch_fast_coef(errorStatus, coefs%coef%solar, coefs%coef%fmv_gas_id, coefs%coef%fmv_coe, coefs%coef%fmv_model_ver, &
4652 coefs%coef%nlayers, coefs%coef%fmv_gas)
4653 if (coefs%coef%fmv_model_ver > 9) THEN
4654 call dispatch_fast_coef(errorStatus, coefs%coef%solar_corr, coefs%coef%fmv_gas_id, coefs%coef%fmv_ncorr, coefs%coef%fmv_model_ver, &
4655 coefs%coef%nlayers, coefs%coef%fmv_gas)
4656 end if
4657 end if
4658
4659 if (coefs%coef%nltecoef) then
4660
4661 call rpn_comm_bcast(coefs%coef%nlte_coef%ncoef, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4662 call rpn_comm_bcast(coefs%coef%nlte_coef%nsol, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4663 call rpn_comm_bcast(coefs%coef%nlte_coef%nsat, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4664 call rpn_comm_bcast(coefs%coef%nlte_coef%nchan, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4665 call rpn_comm_bcast(coefs%coef%nlte_coef%start_chan, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4666
4667 allocate(nlte_chans(coefs%coef%fmv_chn)) ! Index of selected channels in nlte_coefs array in the file
4668
4669 nlte_count = 0 ! Number of NLTE channels being read in
4670 nlte_start = 0 ! Index (in input channel list) of first NLTE channel being read in
4671 do i = 1, coefs%coef%fmv_chn
4672 if (channels(i) >= coefs%coef%nlte_coef%start_chan .and. &
4673 channels(i) < coefs%coef%nlte_coef%start_chan + coefs%coef%nlte_coef%nchan) then
4674 nlte_count = nlte_count + 1
4675 nlte_chans(nlte_count) = channels(i) - coefs%coef%nlte_coef%start_chan + 1
4676 if (nlte_count == 1) nlte_start = i
4677 end if
4678 end do
4679
4680 coefs%coef%nltecoef = ( nlte_count > 0)
4681
4682 nlte_file_nchan = coefs%coef%nlte_coef%nchan
4683
4684 ! Reset NLTE channel variables according to input channel limit
4685 coefs%coef%nlte_coef%start_chan = nlte_start
4686 coefs%coef%nlte_coef%nchan = nlte_count
4687
4688 if (mmpi_myid > 0) then
4689 allocate (coefs%coef%nlte_coef%sec_sat(coefs%coef%nlte_coef%nsat) )
4690 allocate (coefs%coef%nlte_coef%sol_zen_angle(coefs%coef%nlte_coef%nsol) )
4691 allocate (coefs%coef%nlte_coef%sat_zen_angle(coefs%coef%nlte_coef%nsat) )
4692 allocate (coefs%coef%nlte_coef%cos_sol(coefs%coef%nlte_coef%nsol) )
4693 end if
4694 call rpn_comm_bcast(coefs%coef%nlte_coef%sec_sat, coefs%coef%nlte_coef%nsat, 'MPI_REAL8', 0, 'GRID', ierr)
4695 call rpn_comm_bcast(coefs%coef%nlte_coef%sol_zen_angle, coefs%coef%nlte_coef%nsol, 'MPI_REAL8', 0, 'GRID', ierr)
4696 call rpn_comm_bcast(coefs%coef%nlte_coef%sat_zen_angle, coefs%coef%nlte_coef%nsat, 'MPI_REAL8', 0, 'GRID', ierr)
4697 call rpn_comm_bcast(coefs%coef%nlte_coef%cos_sol, coefs%coef%nlte_coef%nsol, 'MPI_REAL8', 0, 'GRID', ierr)
4698
4699 allocate(bigArray(coefs%coef%nlte_coef%ncoef, coefs%coef%nlte_coef%nsat, &
4700 coefs%coef%nlte_coef%nsol, nlte_file_nchan) )
4701
4702 if (mmpi_myid == 0) then
4703 do ichan = 1, nlte_file_nchan
4704 do isol = 1, coefs%coef%nlte_coef%nsol
4705 do isat = 1, coefs%coef%nlte_coef%nsat
4706 do I=1, coefs%coef%nlte_coef%ncoef
4707 bigArray(i, isat, isol, ichan) = coefs%coef%nlte_coef%coef(i,isat,isol,ichan)
4708 end do
4709 end do
4710 end do
4711 end do
4712 deallocate( coefs%coef%nlte_coef%coef )
4713 end if
4714
4715 call rpn_comm_bcast(bigArray, size(bigArray), 'MPI_REAL8', 0, 'GRID', ierr)
4716 allocate(coefs%coef%nlte_coef%coef(coefs%coef%nlte_coef%ncoef, coefs%coef%nlte_coef%nsat, &
4717 coefs%coef%nlte_coef%nsol, nlte_count))
4718 coefs%coef%nlte_coef%coef(:,:,:,:) = bigArray(:,:,:,nlte_chans(1:nlte_count))
4719 deallocate(nlte_chans, bigArray)
4720
4721 end if
4722
4723 if (mmpi_myid==0 .and. associated(coefs%coef%bounds) ) deallocate(coefs%coef%bounds)
4724
4725 !allocate bounds array to store opdep calculation layer limits
4726 !1st dim: upper boundary layer [ub](above which coefs all zeros), lower boundary layer [lb]
4727 !4th dim: thermal layer limits, solar layer limits
4728 allocate(coefs%coef%bounds(2, coefs%coef%fmv_gas, coefs%coef%fmv_chn, 2))
4729 call set_fastcoef_level_bounds(coefs%coef, coefs%coef%thermal, thermal = .true._jplm)
4730 ! if the solar_fast_coefficients section is not present then point the solar coefs to the thermal coefs
4731 if (coefs%coef%solarcoef) then
4732 call set_fastcoef_level_bounds(coefs%coef, coefs%coef%solar, thermal = .false._jplm)
4733 else
4734 coefs%coef%solar => coefs%coef%thermal
4735 coefs%coef%solar_corr => coefs%coef%thermal_corr
4736 coefs%coef%bounds(:,:,:,2) = coefs%coef%bounds(:,:,:,1)
4737 end if
4738
4739 coefs%coef%ff_val_bc = any(coefs%coef%ff_bco(:) /= 0.0d0) .or. any(coefs%coef%ff_bcs(:) /= 1.d0)
4740 coefs%coef%ff_val_gam = any(coefs%coef%ff_gam(:) /= 1.d0)
4741
4742 ! surface water reflectance for visible/near-ir channels
4743 if (any(coefs%coef%ss_val_chn == 2)) then
4744 if ( mmpi_myid==0) deallocate(coefs%coef%refl_visnir_ow, &
4745 coefs%coef%refl_visnir_fw, stat = errorStatus)
4746 allocate(coefs%coef%refl_visnir_ow(coefs%coef%fmv_chn), &
4747 coefs%coef%refl_visnir_fw(coefs%coef%fmv_chn), stat = errorStatus)
4748 call rttov_refl_water_interp(coefs%coef%ff_cwn, coefs%coef%refl_visnir_ow, coefs%coef%refl_visnir_fw)
4749 end if
4750
4751 if (coefs%coef%pmc_shift .and. mmpi_myid > 0) then
4752 allocate(coefs%coef%pmc_ppmc(coefs%coef%fmv_chn), stat = errorStatus)
4753 else
4754 nullify(coefs%coef%pmc_pnominal, coefs%coef%pmc_coef, coefs%coef%pmc_ppmc)
4755 end if
4756
4757 contains
4758
4759 subroutine nullify_gas_coef_pointers(fast_coef)
4760 implicit none
4761
4762 ! Arguments:
4763 type(rttov_fast_coef), intent(inout) :: fast_coef
4764
4765 nullify (fast_coef%mixedgas,&
4766 fast_coef%watervapour, &
4767 fast_coef%ozone, &
4768 fast_coef%wvcont, &
4769 fast_coef%co2, &
4770 fast_coef%n2o, &
4771 fast_coef%co, &
4772 fast_coef%ch4, &
4773 fast_coef%so2)
4774 end subroutine nullify_gas_coef_pointers
4775
4776 subroutine dispatch_fast_coef(err, fast_coef, gas_ids, ncoefs, version, nlayers, ngas)
4777 implicit none
4778
4779 ! Arguments:
4780 integer, intent(out) :: err
4781 type(rttov_fast_coef), pointer, intent(inout) :: fast_coef(:)
4782 integer(jpim), intent(in) :: gas_ids(:)
4783 integer(jpim), intent(in) :: ncoefs(:)
4784 integer(jpim), intent(in) :: version
4785 integer(jpim), intent(in) :: nlayers
4786 integer(jpim), intent(in) :: ngas
4787
4788 ! Locals:
4789 integer(jpim) :: channelIndex, gasIndex, layerIndex, coefIndex
4790 real(8), allocatable :: bigArray(:,:,:,:)
4791 logical :: allocated0
4792
4793 allocate(bigArray(countUniqueChannel,maxval(ncoefs),ngas,nlayers), stat=err )
4794 bigArray(:,:,:,:) = 0.0d0
4795
4796 if (mmpi_myid > 0) then
4797 allocate (fast_coef(countUniqueChannel) )
4798 do channelIndex = 1, countUniqueChannel
4799 allocate(fast_coef(channelIndex)%gasarray(ngas) )
4800 do gasIndex = 1, ngas
4801 allocate (fast_coef(channelIndex)%gasarray(gasIndex)%coef( ncoefs(gasIndex), nlayers) )
4802 end do
4803 end do
4804 end if
4805
4806 do channelIndex = 1, countUniqueChannel
4807 do gasIndex = 1, ngas
4808 call broadcastR82dArray( fast_coef(channelIndex)%gasarray(gasIndex)%coef )
4809 end do
4810 end do
4811
4812 do channelIndex = 1, countUniqueChannel
4813 do gasIndex = 1, ngas
4814 associated0 = associated( fast_coef(channelIndex)%gasarray(gasIndex)%coef )
4815 call rpn_comm_bcast(associated0, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
4816 if (associated0) then
4817 do layerIndex=1, nlayers
4818 do coefIndex=1, ncoefs(gasIndex)
4819 bigArray(channelIndex,coefIndex,gasIndex,layerIndex) = fast_coef(channelIndex)%gasarray(gasIndex)%coef(coefIndex,layerIndex)
4820 end do
4821 end do
4822 end if
4823 end do
4824 end do
4825
4826 do channelIndex = 1, countUniqueChannel
4827 do gasIndex = 1, ngas
4828 associated0 = associated(fast_coef(channelIndex)%gasarray(gasIndex)%coef)
4829 if (associated0) deallocate(fast_coef(channelIndex)%gasarray(gasIndex)%coef)
4830 end do
4831 deallocate(fast_coef(channelIndex)%gasarray)
4832 end do
4833 deallocate(fast_coef)
4834 allocate (fast_coef(coefs%coef%fmv_chn) )
4835 do channelIndex = 1, coefs%coef%fmv_chn
4836 allocate (fast_coef(channelIndex)%gasarray(ngas) )
4837 call nullify_gas_coef_pointers( fast_coef(channelIndex) )
4838 do gasIndex = 1, ngas
4839 if (any( bigArray(indexchan(channelIndex),:,gasIndex,:) /= 0.) ) then
4840 allocate (fast_coef(channelIndex)%gasarray(gasIndex)%coef( ncoefs(gasIndex), nlayers) )
4841 do layerIndex=1, nlayers
4842 do coefIndex=1, ncoefs(gasIndex)
4843 fast_coef(channelIndex)%gasarray(gasIndex)%coef(coefIndex,layerIndex) = bigArray(indexchan(channelIndex),coefIndex,gasIndex,layerIndex)
4844 end do
4845 end do
4846 end if
4847 call set_pointers(fast_coef(channelIndex), gasIndex, gas_ids(gasIndex))
4848 end do
4849 end do
4850
4851 deallocate(bigArray, stat=err )
4852
4853 end subroutine dispatch_fast_coef
4854
4855 end subroutine tvs_rttov_read_coefs
4856
4857
4858 subroutine extractI41dArray(array,oldSize,index)
4859 implicit none
4860
4861 ! Arguments:
4862 integer, pointer, intent(inout) :: array(:)
4863 integer, intent(in) :: oldSize
4864 integer, intent(in) :: index(:)
4865
4866 ! Locals:
4867 integer :: newSize, tmpI41d(oldSize), ierr, trueSize
4868
4869 if (mmpi_myid == 0) then
4870 if (associated(array)) then
4871 trueSize = size(array)
4872 else
4873 trueSize = 0
4874 end if
4875 end if
4876 ierr = 0
4877 call rpn_comm_bcast(trueSize, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4878 if (ierr /= 0) then
4879 write(*,*) 'error 1 in rpn_comm_bcast', ierr, trueSize
4880 call utl_abort('extractI41dArray')
4881 end if
4882 if (trueSize < 1) return
4883
4884 if (trueSize /= oldSize) then
4885 write(*,*) 'extractI41dArray: should not happen ', trueSize, oldSize
4886 end if
4887
4888 newSize = size( index )
4889
4890 if (mmpi_myid > 0) allocate( array(oldSize))
4891 ierr = 0
4892 call rpn_comm_bcast(array, oldSize, 'MPI_INTEGER', 0, 'GRID', ierr)
4893 if (ierr /= 0) then
4894 write(*,*) 'error 2 in rpn_comm_bcast', ierr, array(:)
4895 call utl_abort('extractI41dArray')
4896 end if
4897 tmpI41d = array
4898 deallocate( array )
4899 allocate( array(newSize))
4900 array( : ) = tmpI41d( index(:) )
4901 end subroutine extractI41dArray
4902
4903 subroutine extractR81dArray(array,oldSize,index)
4904 implicit none
4905
4906 ! Arguments:
4907 real(8), pointer, intent(inout) :: array(:)
4908 integer, intent(in) :: oldSize
4909 integer, intent(in) :: index(:)
4910
4911 ! Locals:
4912 integer :: newSize, ierr, trueSize
4913 real(8) :: tmpR81d(oldSize)
4914
4915 if (mmpi_myid == 0) then
4916 if (associated(array)) then
4917 trueSize = size(array)
4918 else
4919 trueSize = 0
4920 end if
4921 end if
4922 ierr = 0
4923 call rpn_comm_bcast(trueSize, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4924 if (ierr /= 0) then
4925 write(*,*) 'error 1 in rpn_comm_bcast', ierr, trueSize
4926 call utl_abort('extractR81dArray')
4927 end if
4928 if (trueSize < 1) return
4929
4930 if (trueSize /= oldSize) then
4931 write(*,*) 'extractR81dArray: should not happen ', trueSize, oldSize
4932 end if
4933
4934 newSize = size( index )
4935
4936 if (mmpi_myid > 0) allocate( array(oldSize))
4937 ierr = 0
4938 call rpn_comm_bcast(array, oldSize, 'MPI_REAL8', 0, 'GRID', ierr)
4939 if (ierr /= 0) then
4940 write(*,*) 'error 2 in rpn_comm_bcast', ierr, array(:)
4941 call utl_abort('extractR81dArray')
4942 end if
4943 tmpR81d = array
4944 deallocate( array )
4945 allocate( array(newSize))
4946 array( : ) = tmpR81d( index(:) )
4947 end subroutine extractR81dArray
4948
4949
4950 subroutine extractR82dArray(array,oldSize1,oldSize2,index)
4951 !second dimension is for channels
4952 implicit none
4953
4954 ! Arguments:
4955 real(8), pointer, intent(inout) :: array(:,:)
4956 integer, intent(in) :: oldSize1
4957 integer, intent(in) :: oldSize2
4958 integer, intent(in) :: index(:)
4959
4960 ! Locals:
4961 integer :: newSize, ierr, trueSize,i
4962 real(8) :: tmpR82d(oldSize1,oldsize2)
4963
4964 if (mmpi_myid == 0) then
4965 if (associated(array)) then
4966 trueSize = size(array)
4967 else
4968 trueSize = 0
4969 end if
4970 end if
4971 ierr = 0
4972 call rpn_comm_bcast(trueSize, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
4973 if (ierr /= 0) then
4974 write(*,*) 'error 1 in rpn_comm_bcast', ierr, trueSize
4975 call utl_abort('extractR82dArray')
4976 end if
4977 if (trueSize < 1) return
4978
4979 if (trueSize /= oldSize1 * oldSize2) then
4980 write(*,*) 'extractR82dArray: should not happen ', trueSize, oldSize1, oldSize2
4981 end if
4982
4983 newSize = size( index )
4984
4985 if (mmpi_myid > 0) allocate( array(oldSize1,oldSize2) )
4986 ierr = 0
4987 call rpn_comm_bcast(array, oldSize1*oldSize2, 'MPI_REAL8', 0, 'GRID', ierr)
4988 if (ierr /= 0) then
4989 write(*,*) 'error 2 in rpn_comm_bcast', ierr, array(:,:)
4990 call utl_abort('extractR82dArray')
4991 end if
4992 tmpR82d = array
4993 deallocate( array )
4994 allocate( array(oldSize1,newSize))
4995 do i=1, newSize
4996 array( :,i) = tmpR82d(:, index(i) )
4997 end do
4998 end subroutine extractR82dArray
4999
5000 subroutine extractR83dArray(array,oldSize1,oldSize2,oldSize3,index)
5001 !second dimension is for channels
5002 implicit none
5003
5004 ! Arguments:
5005 real(8), pointer, intent(inout) :: array(:,:,:)
5006 integer, intent(in) :: oldSize1
5007 integer, intent(in) :: oldSize2
5008 integer, intent(in) :: oldSize3
5009 integer, intent(in) :: index(:)
5010
5011 ! Locals:
5012 integer :: newSize, ierr, trueSize,i
5013 real(8) :: tmpR83d(oldSize1,oldSize2,oldSize3)
5014
5015 if (mmpi_myid == 0) then
5016 if (associated(array)) then
5017 trueSize = size(array)
5018 else
5019 trueSize = 0
5020 end if
5021 end if
5022 ierr = 0
5023 call rpn_comm_bcast(trueSize, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
5024 if (ierr /= 0) then
5025 write(*,*) 'error 1 in rpn_comm_bcast', ierr, trueSize
5026 call utl_abort('extractR83dArray')
5027 end if
5028 if (trueSize < 1) return
5029
5030 if (trueSize /= oldSize1 * oldSize2 * oldSize3) then
5031 write(*,*) 'extractR83dArray: should not happen ', trueSize, oldSize1, oldSize2, oldSize3
5032 end if
5033
5034 newSize = size( index )
5035
5036 if (mmpi_myid > 0) allocate( array(oldSize1,oldSize2, oldSIze3) )
5037 ierr = 0
5038 call rpn_comm_bcast(array, oldSize1*oldSize2*oldSize3, 'MPI_REAL8', 0, 'GRID', ierr)
5039 if (ierr /= 0) then
5040 write(*,*) 'error 2 in rpn_comm_bcast', ierr, array(:,:,:)
5041 call utl_abort('extractR83dArray')
5042 end if
5043 tmpR83d = array
5044 deallocate( array )
5045 allocate( array(oldSize1,newSize,oldSize3))
5046 do i=1, newSize
5047 array( :,i,:) = tmpR83d(:, index(i),: )
5048 end do
5049 end subroutine extractR83dArray
5050
5051 subroutine extractCmplx81dArray(array,oldSize,index)
5052 implicit none
5053
5054 ! Arguments:
5055 complex(kind=8), pointer, intent(inout) :: array(:)
5056 integer, intent(in) :: oldSize
5057 integer, intent(in) :: index(:)
5058
5059 ! Locals:
5060 integer :: newSize, ierr, trueSize
5061 complex(kind=8) :: tmpCx81d(oldSize)
5062
5063 if (mmpi_myid == 0) then
5064 if (associated(array)) then
5065 trueSize = size(array)
5066 else
5067 trueSize = 0
5068 end if
5069 end if
5070 ierr = 0
5071 call rpn_comm_bcast(trueSize, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
5072 if (ierr /= 0) then
5073 write(*,*) 'error 1 in rpn_comm_bcast', ierr, trueSize
5074 call utl_abort('extractCmplx81dArray')
5075 end if
5076 if (trueSize < 1) return
5077
5078 if (trueSize /= oldSize) then
5079 write(*,*) 'extractCmplx81dArray: should not happen ', trueSize, oldSize
5080 end if
5081
5082 newSize = size( index )
5083
5084 if (mmpi_myid > 0) allocate( array(oldSize))
5085 ierr = 0
5086 call rpn_comm_bcast(array, oldSize, 'MPI_COMPLEX8', 0, 'GRID', ierr)
5087 if (ierr /= 0) then
5088 write(*,*) 'error 2 in rpn_comm_bcast', ierr, array(:)
5089 call utl_abort('extractCmplx81dArray')
5090 end if
5091 tmpCx81d = array
5092 deallocate( array )
5093 allocate( array(newSize))
5094 array( : ) = tmpCx81d( index(:) )
5095 end subroutine extractCmplx81dArray
5096
5097
5098 subroutine broadcastR82dArray(array)
5099 implicit none
5100
5101 ! Arguments:
5102 real(kind=8), pointer, intent(inout) :: array(:,:)
5103
5104 ! Locals:
5105 logical :: associated0
5106 integer :: ierr
5107
5108 associated0 = associated(array)
5109 call rpn_comm_bcast(associated0, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
5110 if (ierr /= 0) then
5111 write(*,*) 'error 1 in rpn_comm_bcast', ierr, associated0
5112 call utl_abort('broadcastR82dArray')
5113 end if
5114 ierr = 0
5115 if (associated0) call rpn_comm_bcast(array, size(array) , 'MPI_REAL8', 0, 'GRID', ierr)
5116 if (ierr /= 0) then
5117 write(*,*) 'error 2 in rpn_comm_bcast', ierr, size(array,dim=1), size(array,dim=2)
5118 call utl_abort('broadcastR82dArray')
5119 end if
5120
5121 end subroutine broadcastR82dArray
5122
5123 subroutine broadcastR81dArray(array)
5124 implicit none
5125
5126 ! Arguments:
5127 real(kind=8), pointer, intent(inout) :: array(:)
5128
5129 ! Locals:
5130 logical :: associated0
5131 integer :: ierr
5132
5133 associated0 = associated(array)
5134 call rpn_comm_bcast(associated0, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
5135 if (ierr/=0) then
5136 write(*,*) 'error 1 in rpn_comm_bcast', ierr, associated0
5137 call utl_abort('broadcastR81dArray')
5138 end if
5139 ierr = 0
5140 if (associated0) call rpn_comm_bcast(array, size(array) , 'MPI_REAL8', 0, 'GRID', ierr)
5141 if (ierr/=0) then
5142 write(*,*) 'error 2 in rpn_comm_bcast', ierr, size(array)
5143 call utl_abort('broadcastR81dArray')
5144 end if
5145
5146 end subroutine broadcastR81dArray
5147
5148
5149 subroutine broadcastI41dArray(array)
5150 implicit none
5151
5152 ! Arguments:
5153 integer(kind=4), pointer, intent(inout) :: array(:)
5154
5155 ! Locals:
5156 logical :: associated0
5157 integer :: ierr
5158
5159 associated0 = associated(array)
5160 call rpn_comm_bcast(associated0, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
5161 if (ierr/=0) then
5162 write(*,*) 'error 1 in rpn_comm_bcast', ierr, associated0
5163 call utl_abort('broadcastI41dArray')
5164 end if
5165 ierr = 0
5166 if (associated0) call rpn_comm_bcast(array, size(array) , 'MPI_INTEGER', 0, 'GRID', ierr)
5167 if (ierr/=0) then
5168 write(*,*) 'error 2 in rpn_comm_bcast', ierr, size(array)
5169 call utl_abort('broadcastI41dArray')
5170 end if
5171
5172 end subroutine broadcastI41dArray
5173
5174 !--------------------------------------------------------------------------
5175 ! tvs_printDetailledOmfStatistics
5176 !--------------------------------------------------------------------------
5177 subroutine tvs_printDetailledOmfStatistics(obsSpaceData)
5178 !
5179 !:Purpose: Print channel by channnel O-F statistics fro radiances
5180 !
5181 implicit none
5182
5183 ! Arguments:
5184 type(struct_obs), intent(inout) :: obsSpaceData! obsSpacaData structure
5185
5186 ! Locals:
5187 integer :: sensorIndex, channelIndex, tovsIndex
5188 real(8) zjoch (0:tvs_maxChannelNumber,tvs_maxNumberOfSensors)
5189 real(8) zavgnrm(0:tvs_maxChannelNumber,tvs_maxNumberOfSensors)
5190 real(pre_obsReal) :: zdtb, obsPRM
5191 integer nchanperline, startChannel, endChannel
5192 integer count, incanjo
5193 integer idatyp
5194 integer rttovChannelNumber, bufrChannelNumber
5195 integer inobsch(0:tvs_maxChannelNumber,tvs_maxNumberOfSensors)
5196 integer lcanjo(tvs_maxChannelNumber)
5197 integer :: headerIndex, bodyIndex
5198 real(8) :: sigmaObs
5199
5200 write(*,*) 'tvs_printDetailledOmfStatistics: Starting'
5201
5202 if ( tvs_nobtov == 0) return ! exit if there are not tovs data
5203
5204 ! 1. Computation of (hx - z)/sigma for tovs data only
5205
5206 count = 0
5207 inobsch(:,:) = 0
5208 zjoch (:,:) = 0.0d0
5209 zavgnrm(:,:) = 0.0d0
5210
5211 ! loop over all header indices of the 'TO' family
5212 call obs_set_current_header_list(obsSpaceData,'TO')
5213 HEADER: do
5214 headerIndex = obs_getHeaderIndex(obsSpaceData)
5215 if (headerIndex < 0) exit HEADER
5216
5217 ! 1.1 Extract general information for this observation point
5218 ! ------------------------------------------------------
5219
5220 ! process only radiance data to be assimilated?
5221 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
5222 if ( .not. tvs_isIdBurpTovs(idatyp) ) then
5223 write(*,*) 'tvs_printDetailledOmfStatistics: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
5224 cycle HEADER
5225 end if
5226 tovsIndex = tvs_tovsIndex(headerIndex)
5227 if ( tovsIndex == -1 ) cycle HEADER
5228
5229 sensorIndex = tvs_lsensor(tovsIndex)
5230
5231 ! Set the body list
5232 ! (& start at the beginning of the list)
5233 call obs_set_current_body_list(obsSpaceData, headerIndex)
5234 count = 0
5235 BODY: do
5236 bodyIndex = obs_getBodyIndex(obsSpaceData)
5237 if (bodyIndex < 0) exit BODY
5238
5239 ! Only consider if flagged for assimilation
5240 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) /= obs_assimilated ) cycle BODY
5241
5242 call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
5243 rttovChannelNumber, channelIndex )
5244 bufrChannelNumber = rttovChannelNumber + tvs_channelOffset(sensorIndex)
5245 if ( channelIndex == 0 ) then
5246 write(*,'(A)') ' tvs_printDetailledOmfStatistics: error with channel number'
5247 call utl_abort(' tvs_printDetailledOmfStatistics')
5248 end if
5249
5250 zdtb = obs_bodyElem_r(obsSpaceData,OBS_PRM,bodyIndex) - &
5251 tvs_radiance (tovsIndex) % bt(channelIndex)
5252 if ( tvs_debug ) then
5253 obsPRM = obs_bodyElem_r(obsSpaceData,OBS_PRM,bodyIndex)
5254 write(*,'(a,i4,2f8.2,f6.2)') ' rttovChannelNumber,sim,obs,diff= ', &
5255 rttovChannelNumber, tvs_radiance (tovsIndex) % bt(channelIndex), &
5256 obsPRM, -zdtb
5257 end if
5258
5259 sigmaObs = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
5260
5261 if ( sigmaObs == MPC_missingValue_R8) cycle body
5262
5263 count = count + 1
5264 inobsch(bufrChannelNumber,sensorIndex) = inobsch(bufrChannelNumber,sensorIndex) + 1
5265 zjoch(bufrChannelNumber,sensorIndex) = &
5266 zjoch(bufrChannelNumber,sensorIndex) &
5267 + zdtb * zdtb / (sigmaObs * sigmaObs)
5268 zavgnrm(bufrChannelNumber,sensorIndex) = &
5269 zavgnrm(bufrChannelNumber,sensorIndex) - &
5270 zdtb / sigmaObs
5271 end do BODY
5272
5273 end do HEADER
5274
5275 ! 2. Close up, print summary
5276 ! . -----------------------
5277
5278
5279 ! printout of mean jo and normalized average for each sensor.
5280
5281 nchanperline = 18
5282 if ( count > 0 ) then
5283 write(*,*)
5284 write(*,*)
5285 write(*,'(10x,A)') '- tvs_printDetailledOmfStatistics: computing jo and residuals to tovs observations'
5286
5287 do sensorIndex = 1, tvs_nsensors
5288 inobsch(0,sensorIndex) = sum ( inobsch(1:,sensorIndex) )
5289 zjoch(0,sensorIndex) = sum( zjoch(1:,sensorIndex) )
5290 zavgnrm(0,sensorIndex) = sum( zavgnrm(1:,sensorIndex) )
5291 end do
5292
5293 do sensorIndex = 1, tvs_nsensors
5294 incanjo = 0
5295 do channelIndex = 0, tvs_maxChannelNumber
5296 if ( inobsch(channelIndex, sensorIndex) /= 0 ) then
5297 incanjo = incanjo + 1
5298 lcanjo(incanjo) = channelIndex
5299 end if
5300 end do
5301 if ( incanjo /= 0 ) then
5302 write(*,"(/1x,'sensor #',i2,'. platform: ',a, 'instrument: ',a)") &
5303 sensorIndex, tvs_satelliteName(sensorIndex), tvs_instrumentName(sensorIndex)
5304 do startChannel = 1, incanjo, nchanperline
5305 endChannel = min(startChannel + nchanperline - 1 , incanjo)
5306 if ( startChannel == 1 ) then
5307 write(*,"(1x,'channel',t13,' all',17i6)") (lcanjo(channelIndex), channelIndex=startChannel+1, endChannel)
5308 else
5309 write(*,"(1x,'channel',t13,18i6)") (lcanjo(channelIndex), channelIndex=startChannel, endChannel)
5310 end if
5311 write(*,"(1x,'no. obs.',t13,18i6)") (inobsch(lcanjo(channelIndex),sensorIndex), channelIndex=startChannel, endChannel)
5312 write(*,"(1x,'mean jo',t13,18f6.2)") &
5313 (zjoch(lcanjo(channelIndex),sensorIndex)/max(1,inobsch(lcanjo(channelIndex),sensorIndex)), channelIndex=startChannel,endChannel)
5314 write(*,"(1x,'norm. bias',t13,18f6.2,/)") &
5315 (zavgnrm(lcanjo(channelIndex),sensorIndex)/max(1,inobsch(lcanjo(channelIndex), sensorIndex)) , channelIndex=startChannel, endChannel)
5316 end do
5317 end if
5318 end do
5319 end if
5320
5321 end subroutine tvs_printDetailledOmfStatistics
5322
5323
5324 !--------------------------------------------------------------------------
5325 ! tvs_getLocalChannelIndexFromChannelNumber
5326 !--------------------------------------------------------------------------
5327 subroutine tvs_getLocalChannelIndexFromChannelNumber(idsat,channelIndex_out,channelNumber_in)
5328 !
5329 !:Purpose: to get local channel index from channel number
5330 !
5331 implicit none
5332
5333 ! Arguments:
5334 integer, intent(in) :: idsat ! Satellite index
5335 integer, intent(out) :: channelIndex_out ! Channel index
5336 integer, intent(in) :: channelNumber_in ! Channel number
5337
5338 ! Locals:
5339 logical, save :: first =.true.
5340 integer :: channelNumber, sensorIndex, channelIndex
5341 integer, allocatable, save :: savedChannelIndexes(:,:)
5342
5343 if (first) then
5344 allocate( savedChannelIndexes(tvs_nsensors, tvs_maxChannelNumber ) )
5345 savedChannelIndexes(:,:) = -1
5346 do sensorIndex = 1, tvs_nsensors
5347 channels:do channelNumber = 1, tvs_maxChannelNumber
5348 indexes: do channelIndex =1, tvs_nchan(sensorIndex)
5349 if ( channelNumber == tvs_ichan(channelIndex,sensorIndex) ) then
5350 savedChannelIndexes(sensorIndex,channelNumber) = channelIndex
5351 exit indexes
5352 end if
5353 end do indexes
5354 end do channels
5355 end do
5356 first = .false.
5357 end if
5358
5359 channelIndex_out = savedChannelIndexes(idsat,channelNumber_in)
5360
5361 if (channelIndex_out == -1) then
5362 write(*,*) 'channel number requested = ', channelNumber_in
5363 write(*,*) 'idsat = ', idsat
5364 write(*,*) 'tvs_getLocalChannelIndexFromChannelNumber: warning channel not found'
5365 end if
5366
5367 end subroutine tvs_getLocalChannelIndexFromChannelNumber
5368
5369
5370 subroutine updateCloudInTovsProfile(sensorTovsIndexes, nlv_T, mode, beSilent)
5371 !
5372 !:Purpose: Modify the cloud in tvs_profiles_nl structure of rttov.
5373 !
5374 implicit none
5375
5376 ! Arguments:
5377 integer, intent(in) :: sensorTovsIndexes(:)
5378 integer, intent(in) :: nlv_T
5379 character(*), intent(in) :: mode ! save or restore
5380 logical, intent(in) :: beSilent ! flag to control verbosity
5381
5382 ! Locals:
5383 integer :: profileIndex, profileCount
5384 real(8), allocatable, save :: cloudProfileToStore(:,:)
5385
5386 if ( .not. beSilent ) write(*,*) 'updateCloudInTovsProfile: Starting'
5387 if ( .not. beSilent ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
5388
5389 profileCount = size(sensorTovsIndexes)
5390
5391 if ( trim(mode) == 'save' ) then
5392 if (allocated(cloudProfileToStore)) deallocate(cloudProfileToStore)
5393 allocate(cloudProfileToStore(nlv_T,profileCount))
5394
5395 do profileIndex = 1, profileCount
5396 cloudProfileToStore(:,profileIndex) = tvs_profiles_nl(sensorTovsIndexes(profileIndex)) % clw(:)
5397 tvs_profiles_nl(sensorTovsIndexes(profileIndex)) % clw(:) = qlim_getMinValueCloud('LWCR')
5398 end do
5399
5400 else if ( trim(mode) == 'restore' ) then
5401 do profileIndex = 1, profileCount
5402 tvs_profiles_nl(sensorTovsIndexes(profileIndex)) % clw(:) = cloudProfileToStore(:,profileIndex)
5403 end do
5404
5405 deallocate(cloudProfileToStore)
5406
5407 else
5408 call utl_abort('updateCloudInTovsProfile: mode should be either "save" or "restore"')
5409
5410 end if
5411
5412 end subroutine updateCloudInTovsProfile
5413
5414
5415 subroutine updateCloudInTovsCloudProfile(sensorTovsIndexes, nlv_T, mode, beSilent)
5416 !
5417 !:Purpose: Modify the cloud in tvs_cld_profiles_nl structure of rttovScatt.
5418 !
5419 implicit none
5420
5421 ! Arguments:
5422 integer, intent(in) :: sensorTovsIndexes(:)
5423 integer, intent(in) :: nlv_T
5424 character(*), intent(in) :: mode ! save or restore
5425 logical, intent(in) :: beSilent ! flag to control verbosity
5426
5427 ! Locals:
5428 integer :: profileIndex, profileCount
5429 real(8), allocatable, save :: rainFluxProfileToStore(:,:)
5430 real(8), allocatable, save :: snowFluxProfileToStore(:,:)
5431 real(8), allocatable, save :: clwProfileToStore(:,:)
5432 real(8), allocatable, save :: ciwProfileToStore(:,:)
5433 real(8), allocatable, save :: cloudFractionProfileToStore(:,:)
5434
5435 if ( .not. beSilent ) write(*,*) 'updateCloudInTovsCloudProfile: Starting'
5436 if ( .not. beSilent ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
5437
5438 profileCount = size(sensorTovsIndexes)
5439
5440 if ( trim(mode) == 'save' ) then
5441 if (allocated(rainFluxProfileToStore)) deallocate(rainFluxProfileToStore)
5442 if (allocated(snowFluxProfileToStore)) deallocate(snowFluxProfileToStore)
5443 if (allocated(clwProfileToStore)) deallocate(clwProfileToStore)
5444 if (allocated(ciwProfileToStore)) deallocate(ciwProfileToStore)
5445 if (allocated(cloudFractionProfileToStore)) deallocate(cloudFractionProfileToStore)
5446 allocate(rainFluxProfileToStore(nlv_T,profileCount))
5447 allocate(snowFluxProfileToStore(nlv_T,profileCount))
5448 allocate(clwProfileToStore(nlv_T,profileCount))
5449 allocate(ciwProfileToStore(nlv_T,profileCount))
5450 allocate(cloudFractionProfileToStore(nlv_T,profileCount))
5451
5452 do profileIndex = 1, profileCount
5453 rainFluxProfileToStore(:,profileIndex) = tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,1)
5454 snowFluxProfileToStore(:,profileIndex) = tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,2)
5455 clwProfileToStore(:,profileIndex) = tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,4)
5456 ciwProfileToStore(:,profileIndex) = tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,5)
5457 cloudFractionProfileToStore(:,profileIndex) = tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro_frac(:,1)
5458 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,1) = qlim_getMinValueCloud('RF')
5459 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,2) = qlim_getMinValueCloud('SF')
5460 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,4) = qlim_getMinValueCloud('LWCR')
5461 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,5) = qlim_getMinValueCloud('IWCR')
5462 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro_frac(:,1) = qlim_getMinValueCloud('CLDR')
5463 end do
5464
5465 else if ( trim(mode) == 'restore' ) then
5466 do profileIndex = 1, profileCount
5467 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,1) = rainFluxProfileToStore(:,profileIndex)
5468 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,2) = snowFluxProfileToStore(:,profileIndex)
5469 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,4) = clwProfileToStore(:,profileIndex)
5470 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro(:,5) = ciwProfileToStore(:,profileIndex)
5471 tvs_cld_profiles_nl(sensorTovsIndexes(profileIndex)) % hydro_frac(:,1) = cloudFractionProfileToStore(:,profileIndex)
5472 end do
5473
5474 deallocate(rainFluxProfileToStore)
5475 deallocate(snowFluxProfileToStore)
5476 deallocate(clwProfileToStore)
5477 deallocate(ciwProfileToStore)
5478 deallocate(cloudFractionProfileToStore)
5479
5480 else
5481 call utl_abort('updateCloudInTovsCloudProfile: mode should be either "save" or "restore"')
5482
5483 end if
5484
5485 end subroutine updateCloudInTovsCloudProfile
5486
5487 !--------------------------------------------------------------------------
5488 ! tvs_getChannelNumIndexFromPPP
5489 !--------------------------------------------------------------------------
5490 subroutine tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
5491 channelNumber, channelIndex )
5492 !
5493 !:Purpose: Get channel number/index from obs_ppp for TO observations.
5494 !
5495 implicit none
5496
5497 ! Arguments:
5498 type(struct_obs), intent(in) :: obsSpaceData
5499 integer, intent(in) :: headerIndex
5500 integer, intent(in) :: bodyIndex
5501 integer, intent(out) :: channelNumber
5502 integer, intent(out) :: channelIndex
5503
5504 ! Locals:
5505 integer :: tovsIndex, sensorIndex
5506
5507 tovsIndex = tvs_tovsIndex(headerIndex)
5508 sensorIndex = tvs_lsensor(tovsIndex)
5509
5510 channelNumber = nint(obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex))
5511 channelNumber = max( 0 , min( channelNumber , tvs_maxChannelNumber + 1))
5512 channelNumber = channelNumber - tvs_channelOffset(sensorIndex)
5513 channelIndex = utl_findloc(tvs_ichan(:,sensorIndex),channelNumber)
5514
5515 end subroutine tvs_getChannelNumIndexFromPPP
5516
5517end module tovsNL_mod
