RTTOV-WRF gives poor brightness temp with add_clouds = true

[braczka]

Describe the bug

This provides more documentation for the obs_def_rttov13.f90 bug, first reported by Joseph Chan through
the dart email list when performing a perfect model experiment using WRF when generating synthetic brightness temperature observations through the RTTOV radiance model.

Case folder for the bug is provided here: /glade/work/bmraczka/DART/models/wrf/RTTOV_bug

1. Build DART: Make sure mkmf.template is modified to account for RTTOV libraries. Used mkmf.template.rttov.intel as template then modify:

RTTOV=/glade/p/cisl/dares/rttov13_ifort
RTTOV = 13

cd /glade/work/bmraczka/DART/models/wrf/work
modify template input.nml to include obs_def_rttov13_mod.f90 as obs_type_file in &preprocess_nml, and &obs_kind_nml must include HIMAWARI_9_AHI_RADIANCE

Make sure obs_def_rttov13_mod.f90 also includes channel bug fix as well.

Include specific changes in obs_def_rttov13.f90 for this perfect model experiment:

Line 243: ! HIMAWARI_8_AHI_RADIANCE, QTY_BRIGHTNESS_TEMPERATURE
Line 244: ! HIMAWARI_9_AHI_RADIANCE, QTY_BRIGHTNESS_TEMPERATURE

Execute quickbuild.sh in work folder

2. Run the perfect model experiment example :

cd /glade/work/bmraczka/DART/models/wrf/RTTOV_bug

Joseph Chan provided a wrf input file (wrfinput_d01) which has been manipulated to generate a simple cloud
cuboid, using create_artificial_clouds.py. Also obs_seq.in observation template is provided.

See modifications in input.nml for &perfect_model_obs_nml and &model_mod_nml;
this will include add_clouds = .true. within &obs_def_rttov_nml

Generate a link to executable perfect_model_obs in current directory

qsub perfect_model_obs.csh

What was the expected outcome?

Expected outcome was for synthetic obs to have cooler brightness temperatures limited to
the portion of the domain covered by the cuboid cloud as shown
here using plot_obs_seq_IR_data.py:

What actually happened?

Brightness temperatures characteristic of clouds, extend beyond domain of cloud cuboid in a horizontal striping
pattern as shown below:

Error Message

No error messages

Which model(s) are you working with?

WRF3.9.1 and RTTOV13

Version of DART

v10.8.3

Build information

Cheyenne, using intel compiler, see RTTOV build settings above

Have you modified the DART code?

A fix, as first suggested by Joseph Chan was to 'zero out' the cloud water assignment within
obs_def_rttov13.f90 by adding the first line below before entering cloud identification section:

..
..

           runtime % profiles(imem) % cloud(:,:) = 0.0_jprb
         if (.not. is_cumulus) then
            ! stratus
            if (surftype == 0) then
               ! 1, Stratus Continental, STCO (kg/kg), convert levels to layers
               runtime % profiles(imem) % cloud(1,:) = 0.5_jprb*(totalwater(ly1idx) + &
                                                               totalwater(ly2idx))

..
..

Further investigation suggests that the cloud water profile (runtime % profiles(imem) % cloud (:,:)) is initialized properly through rttov_alloc_direct.F90 (/glade/p/cisl/dares/rttov13_ifort/src/main/rttov_alloc_direct.F90). Thus at the beginning of runtime the cloud (:,:) variable is allocated successfully for all cloud types (dimension 1; STCO, STMA, CUCC, CUMA etc.) and vertical levels (dimension 2, 1-42). However, as the code progress through the wrf domain, the bug occurs when the cloud type changes from one vertical profile to the next. In these cases, when the cloud type changes between calls, the current cloud type is assigned the proper water values for the current profile, however, the previous cloud type from the previous profile is retained, and is assigned to the current call -- thus the current profile is assigned both the current cloud water values, but erroneously assigns the previous profile cloud type water values. This erroneous assignment explains why cloud-like brightness temperatures can extend beyond the domain of the cuboid cloud itself. This issue arises, in part, because the code assumes there can only be 1 cloud type per profile, thus the water values are assigned to only 1 cloud type within each profile respectively.

It could be this bug is related to the specific intel compiler environment, as other users have not reported this error. Additional testing should be performed using other compilers (e.g. gfortran, cray). Regardless the zeroing out fix as described above has been demonstrated to work for the intel compiler.

Another note unrelated to the cloud brightness temperature bug is that the wind fetch over the ocean (atmos % wfetch) is never allocated properly within obs_def_rttov13.f90. Thus, the fix is to manually assign a wfetch value to the runtime profile as:

..
..

   if (allocated(atmos % wfetch)) then
      ! Wind fetch over the ocean (m)
      runtime % profiles(imem) % s2m % wfetc = atmos % wfetch(imem)
   else
      runtime % profiles(imem) % s2m % wfetc = 100000.
   end if

..
..

This fix is not ideal and should be revisited.

[hkershaw-brown]

wfetc from rttov13 docs

Real wfetc Wind fetch (m) (length of water over which the wind has blown, typical
value 100000m for open ocean) – only used by sea surface solar BRDF
model.

without hardcoding 100000 runtime % profiles(imem) % s2m % wfetc = 100000.

with addsolar = .false.

 2023/12/12  14:56:56  fatal error in module rttov_check_profiles.F90:0289
     invalid wfetc (wind fetch) (profile number =        1)
 2023/12/12  14:56:56  fatal error in module rttov_direct.F90:0361
     
 ERROR FROM:
  source : obs_def_rttov_mod.f90
  routine: do_forward_model
  message:  RTTOV direct error for platform/satellite/sensor id:          31    
        9          56
  message: ... name: HIMAWARI_9_AHI

with
addsolar = .true.
use_wfetch = .true.

 PE 0: perfect_model_obs: Ready to evaluate up to   14784 observations
 ERROR FROM:
  source : wrf/model_mod.f90
  routine: model_interpolate
  message:  cannot interpolate QTY_WIND_FETCH with the current WRF arrays in sta
 te vector

with
addsolar = .false.
use_wfetch = .false.

runs to completion, answers different to addsolar = .true. + hardcoded wfetc