ORAC-CC / orac

Optimal Retrieval of Aerosol and Cloud
GNU General Public License v3.0
30 stars 20 forks source link

Int_CTP is not called on some occasions #54

Open simonrp84 opened 3 years ago

simonrp84 commented 3 years ago

In the main processor, Int_CTP adjusts the RTM temperature profile by extrapolating the lapse rate above the tropopause, which enables us to put overshoots at approximately the right height and also prevents lower altitude clouds from being placed into the stratosphere.

However, if the First Guess is set via the driver file: Ctrl%FG(IPc, :) = SelmCtrl;SelmCtrl;SelmCtrl for example Then the Int_CTP code is not called, as X_MDAD logic doesn't trigger (it appears the index in this function is always 4/IFr when the FG is set via the driver).

This should be changed, as in most cases we want to interpolate the tropopause temp.

I have a fudged solution in my X_MDAD code:

...
   status = 0

   ! Find 11 and 12 micron indices
   Y_Id = Ctrl%Ind%Y_Id(SPixel%spixel_y_to_ctrl_y_index(1:SPixel%Ind%Ny))
   i_spixel_06 = find_in_array(Y_Id, Ctrl%Ind%Y_Id_legacy(I_legacy_0_6x))
   i_spixel_11 = find_in_array(Y_Id, Ctrl%Ind%Y_Id_legacy(I_legacy_11_x))

   i_spixel_06_solar = find_in_array(SPixel%Ind%YSolar(1:SPixel%Ind%NSolar), &
                                     i_spixel_06)
   call Int_CTP(SPixel, Ctrl, SPixel%Ym(i_spixel_11), X_dump, status)

   ! Parameters supported are Tau, Pc and f.
   select case (index)

   case (ITau, ITau2) ! Cloud optical depth, Tau

      if (SPixel%Illum(1) == IDay .and. i_spixel_06_solar > 0) then
         ! Calculate overcast reflectance (assuming fully cloudy pixel).
         ! Uses channel nearest 0.67 microns, index Ctrl%Ind%MDAD_SW.
         Ref_o = SPixel%Ym(i_spixel_06) - SPixel%Surface%Rs(i_spixel_06_solar)

         ! Convert albedo (range 0 - 1) into index (range 1 to 10)
         ! Use the first SEC_o value, assuming that all values are quite close
         ! and we only need an approximation for first guess setting.
         iFGOP = int( ( Ref_o * SPixel%Geom%SEC_o(1) * 10.0 ) + 1.5 )

         if (iFGOP > 11) then
            iFGOP = 11
         else if (iFGOP < 1) then
            iFGOP = 1
         end if

         X = FGOP(iFGOP)
         if (present(Err)) Err = MDADErrTau
      else ! Can't calculate Tau unless it's daylight
         write(*,*)'ERROR: X_MDAD(): Cant calculate Tau unless its daylight'
         status = XMDADMeth
      end if

   case (IPc, IPc2) ! Cloud pressure, Pc

      if (i_spixel_11 > 0) then
         if (SPixel%Ym(i_spixel_11) /= MissingXn) then
            ! Interpolate for the BT to the rad. profile to get Pc FG/AP
            !call Int_CTP(SPixel, Ctrl, SPixel%Ym(i_spixel_11), X, status)
            X = X_dump
            if (present(Err)) Err = MDADErrPc
         else
            ! Invalid data available
            status = XMDADMeth
            write(*,*) 'WARNING: X_MDAD(): Invalid thermal data'
         end if
      else ! Can't calculate Pc if required LW channels not selected
         status = XMDADMeth
!         write(*,*) 'WARNING: X_MDAD(): Cant calculate Pc if required LW ' // &
!                    'channels not selected'
      end if

   case (IFr) ! Cloud fraction, f
      ! Always overcast
      X = 1
      if (present(Err)) Err = MDADErrF
   end select
...

i.e: I moved the Int_CTP call outside the logic and stored the result in an X_dump variable, which - if required - is then put into X in the place that Int_CTP used to be called. This is probably not the most elegant solution but it works. We should fix this properly.