#include "rundeck_opts.h"

      MODULE RAD_COM 36,3
!@sum  RAD_COM Model radiation arrays and parameters
!@auth Original Development Team
!@ver  1.0
      USE MODEL_COM, only : im,jm,lm,lm_req
      USE DOMAIN_DECOMP, only : grid
      USE RADPAR, only : S0,ITRMAX
!@var S0 solar 'constant' needs to be saved between calls to radiation
      IMPLICIT NONE
      SAVE

C**** DEFAULT ORBITAL PARAMETERS FOR EARTH
C**** Note PMIP runs had specified values that do not necesarily
C**** coincide with those used as the default, or the output of ORBPAR.
C****                    OMEGT          OBLIQ        ECCEN
C**** DEFAULT (2000 AD): 282.9          23.44        0.0167
C**** PMIP CONTROL:      282.04         23.446       0.016724
C**** PMIP 6kyr BP:      180.87         24.105       0.018682
C**** PMIP LGM (21k):    294.42         22.949       0.018994
!@param OMEGT_def precession angle (degrees from vernal equinox)
      real*8, parameter :: omegt_def = 282.9d0
!@param OBLIQ_def obliquity angle  (degrees)
      real*8, parameter :: obliq_def = 23.44d0
!@param ECCN_def eccentricity
      real*8, parameter :: eccn_def  = .0167d0
!@var OMEGT,OBLIQ,ECCN actual orbital parameters used
      real*8 OMEGT,OBLIQ,ECCN

C**** Database parameters to control orbital parameter calculation
C**** Note: setting calc_orb_par with paleo_orb_yr=2000 does not produce
C**** exactly the same as the default values.
!@dbparam calc_orb_par = 1 to calc orbital parameters
      integer :: calc_orb_par = 0
!@dbparam paleo_orb_yr is paleo year (BP) for orbital calc
      real*8 :: paleo_orb_yr = -50.  ! (i.e. 2000AD)

!@var dimrad_sv dimension sum of input fields saved for radia_only runs
      INTEGER, PARAMETER :: dimrad_sv=IM*JM*(7*LM+3*LM_REQ+24)
!@var RQT Radiative equilibrium temperatures above model top
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: RQT
!@var Tchg Total temperature change in adjusted forcing runs
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: Tchg
!@var SRHR(0) Solar   raditive net flux into the ground          (W/m^2)
!@var TRHR(0) Thermal raditive net flux into ground(W/O -StB*T^4)(W/m^2)
!@*   Note: -StB*T^4 is MISSING, since T may vary a lot betw. rad. calls
!@var SRHR(1->LM) Solar   raditive heating rate (W/m^2)  (short wave)
!@var TRHR(1->LM) Thermal raditive heating rate (W/m^2)  (long wave)
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: SRHR,TRHR
!@var TRSURF upward thermal radiation at the surface from rad step W/m2
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: TRSURF
!@var FSF Solar Forcing over each type (W/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: FSF
!@var FSRDIR Direct beam solar incident at surface (W/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: FSRDIR
!@var SRVISSURF Incident solar direct+diffuse visible at surface (W/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: SRVISSURF
!@var SRDN Total incident solar at surface (W/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: SRDN  ! saved in rsf
!@var srnflb_save  Net solar radiation (W/m^2)
!@var trnflb_save  Net thermal radiation (W/m^2)
      REAL*8,ALLOCATABLE,DIMENSION(:,:,:) :: srnflb_save,trnflb_save
!@var ttausv_save  Tracer optical thickness
!@var ttausv_cs_save  Tracer optical thickness clear sky
      REAL*8,ALLOCATABLE,DIMENSION(:,:,:,:) :: ttausv_save,
     &     ttausv_cs_save

!@var CFRAC Total cloud fraction as seen be radiation
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: CFRAC ! saved in rsf
!@var RCLD Total cloud optical depth as seen be radiation
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: RCLD ! saved in rsf
!@var O3_tracer_save 3D ozone saved elsewhere for use in radiation
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: O3_tracer_save!saved rsf
!@var rad_to_chem save 3D quantities from radiation code for use in
!@+   chemistry (or rest of model). 1=Ozone, 2=aerosol ext, 3=N2O, 4=CH4,
!@+   5=CFC11+CFC12
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: rad_to_chem !saved in rsf
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: rad_to_file
!@var KLIQ Flag indicating dry(0)/wet(1) atmosphere (memory feature)
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:) :: KLIQ ! saved in rsf
!@dbparam Ikliq 0,1,-1 initialize kliq as dry,equil,current model state
      INTEGER :: Ikliq = -1  !  get kliq-array from restart file
!@dbparam RHfix const.rel.humidity passed to radiation for aeros. tests
      REAL*8 :: RHfix = -1.  !  pass the current model rel.humidity
!@dbparam dalbsnX global coeff for snow alb change by black carbon depos
      REAL*8 ::  dalbsnX = 0.
!@dbparam albsn_yr year of blk carb depos used for snow alb. reduction
      INTEGER ::  albsn_yr = 1951

!     variables related to aerosol indirect effects:
!     (CDNC=cloud droplet number concentration)
!@dbparam CC_CDNCx scaling factor relating cld cvr change and CDNC change
      REAL*8 :: CC_CDNCX = .0000d0  ! .0036d0
!@dbparam OC_CDNCx scaling factor relating cld opt depth and CDNC change
      REAL*8 :: OD_CDNCX = .0000d0  ! .007d0
!@var pcdnc,vcdnc pressure,vertical profile for cld.cvr change
      real*8, parameter, dimension(7) ::
     * pcdnc=(/984.d0, 964.d0, 934.d0, 884.d0, 810.d0, 710.d0, 550.d0/)
     *,vcdnc=(/ .35d0,  .20d0,  .10d0,  .17d0,  .10d0,  .08d0,   0.d0/)
!@var cdncl = vcdnc interpolated to current vertical resolution
      real*8 cdncl(LM)

!@var COSZ1 Mean Solar Zenith angle for curr. physics(not rad) time step
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: COSZ1
!@dbparam S0X solar constant multiplication factor
      REAL*8 :: S0X = 1.
!@dbparam S0_yr,S0_day obs.date of solar constant (if 0: time var)
      INTEGER :: S0_yr = 1951 , S0_day = 182
!@dbparam CO2X,... scaling factors for CO2 N2O CH4 CFC11 CFC12 XGHG
      REAL*8 :: CO2X=1.,N2OX=1.,CH4X=1., CFC11X=1.,CFC12X=1.,XGHGX=1.
!@dbparam GHG_yr,GHG_day obs.date of well-mixed GHgases (if 0: time var)
      INTEGER :: GHG_yr = 1951 , GHG_day = 182
!@dbparam Volc_yr,Volc_day obs.date of Volc.Aerosols (if 0: time var)
      INTEGER :: Volc_yr = 1951 , Volc_day = 182
!@dbparam Aero_yr obs.year of troposph.Aerosols (if 0: time var)
      INTEGER :: Aero_yr = 1951    ! always use annual cycle
!@dbparam O3_yr obs.year of Ozone (if 0: time var, -yyyy: use decyyyy,
!@+   janyyyy to find jan1-15,dec16-31,not dec of yyyy-1, jan of yyyy+1)
      INTEGER :: O3_yr = 1951      ! always use annual cycle
!@dbparam crops_yr obs.year of crops (if 0: time var, -1: default)
      INTEGER :: crops_yr = -1
!@dbparam H2OstratX strat_water_vapor, cloud, Ozone scaling factor
      REAL*8 :: H2OstratX = 1. , cldX = 1. , O3X = 1.
!@dbparam H2ObyCH4 if not 0: add CH4 produced H2O into layers 1->LM
      REAL*8 :: H2ObyCH4 = 1.
!@var dH2O  zonal H2O-prod.rate in kg/m^2/ppm_CH4/second in layer L
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: dH2O
!@var RSDIST,SIND,COSD orbit related variables computed once a day
      REAL*8 :: RSDIST,SIND,COSD
!@var depoBC,depoBC_1990 observed black carbon deposition (curr,1990)
      REAL*8, DIMENSION(72,46) :: depoBC,depoBC_1990  ! rad.grid
!@var ALB is SRNFLB(1)/(SRDFLB(1)+1.D-20),PLAVIS,PLANIR,ALBVIS,ALBNIR,
!@+       SRRVIS,SRRNIR,SRAVIS,SRANIR (see RADIATION)
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: ALB

!@var SALB (1.-broadband surface albedo) - saved in rsf
      REAL*8, POINTER, DIMENSION(:,:) :: SALB   ! = ALB(:,:,1)
!      EQUIVALENCE (SALB,ALB)

#ifdef ALTER_RADF_BY_LAT
!@var FULGAS_lat multiplicative factors for altering FULGAS by latitude
!@+ for non-transient runs. (greenhouse gas regional forcing)
      real*8, dimension(13,46):: FULGAS_lat !rad not model grid, 13 gasses
!@var FS8OPX_lat multiplicative factors for altering FS8OPX by latitude
!@+ for non-transient runs. (aerosol regional forcing) SOLAR
!@var FT8OPX_lat multiplicative factors for altering FT8OPX by latitude
!@+ for non-transient runs. (aerosol regional forcing) THERMAL
      real*8, dimension(8,46):: FS8OPX_lat,FT8OPX_lat !rad not model grid
                                                !8 groups of aerosols
#endif
!@dbparam rad_interact_tr =1 for radiatively active tracers (default=0)
      INTEGER :: rad_interact_tr = 0
C**** the radiative forcing level for instantaneous forcing calcs is set
C**** using the rad_forc_lev parameter.
!@dbparam rad_forc_lev = 0 for TOA, 1 for LTROPO (default=0)
      INTEGER :: rad_forc_lev = 0

!@dbparam cloud_rad_forc = 1 for calculation of cloud radiative forcing
      INTEGER :: cloud_rad_forc = 0

!@var co2ppm Current CO2 level as seen by radiation
      REAL*8 :: co2ppm = 280.    ! set a resaonable default value

C**** Local variables initialised in init_RAD
!@var PLB0,QL0 global parts of local arrays (to avoid OMP-copyin)
      REAL*8, DIMENSION(LM_REQ)       :: PLB0,SHL0
!@var SINJ,COSJ sines and cosines for zenith angle calculation
      REAL*8, ALLOCATABLE, DIMENSION(:) :: SINJ,COSJ

!@var NTRIX Indexing array for optional aerosol interaction
      INTEGER, DIMENSION(ITRMAX) :: NTRIX = 0
!@var WTTR weighting array for optional aerosol interaction
      REAL*8, DIMENSION(ITRMAX) :: WTTR = 1.
!@var nrad_clay index of clay in arrays for optional aerosol interaction
      INTEGER :: nrad_clay

      END MODULE RAD_COM


      SUBROUTINE ALLOC_RAD_COM(grid) 1,6
!@sum  To allocate arrays who sizes now need to be determined at
!@+    run-time
!@auth Rodger Abel
!@ver  1.0

      USE DOMAIN_DECOMP, ONLY : DIST_GRID
      USE DOMAIN_DECOMP, ONLY : GET
      USE MODEL_COM, ONLY : IM, JM, LM, LM_REQ
#ifdef TRACERS_ON
      USE tracer_com,ONLY : Ntm
#endif
      USE RAD_COM, ONLY : RQT,Tchg,SRHR,TRHR,FSF,FSRDIR,SRVISSURF,TRSURF
     *     ,SRDN, CFRAC, RCLD, O3_tracer_save,rad_to_chem,rad_to_file
     *     ,KLIQ, COSZ1, dH2O, ALB, SALB, SINJ, COSJ,srnflb_save,
     *     trnflb_save, ttausv_save, ttausv_cs_save

      IMPLICIT NONE
      TYPE (DIST_GRID), INTENT(IN) :: grid

      INTEGER :: J_0H, J_1H
      INTEGER :: IER

      CALL GET(grid, J_STRT_HALO=J_0H, J_STOP_HALO=J_1H)

      ALLOCATE( RQT(LM_REQ, IM, J_0H:J_1H),
     *     Tchg(LM+LM_REQ, IM, J_0H:J_1H),
     *     SRHR(0:LM, IM, J_0H:J_1H),
     *     TRHR(0:LM, IM, J_0H:J_1H),
     *     TRSURF(1:4, IM, J_0H:J_1H),
     *     FSF(4,IM, J_0H:J_1H),
     *     FSRDIR(IM, J_0H:J_1H),
     *     SRVISSURF(IM, J_0H:J_1H),
     *     SRDN(IM, J_0H:J_1H),
     *     CFRAC(IM, J_0H:J_1H),
     *     RCLD(LM, IM, J_0H:J_1H),
     *     O3_tracer_save(LM, IM, J_0H:J_1H),
     *     rad_to_chem(LM, IM, J_0H:J_1H,5),
     *     rad_to_file(LM, IM, J_0H:J_1H,5),
     *     KLIQ(LM,4, IM, J_0H:J_1H),
     *     COSZ1(IM, J_0H:J_1H),
     *     dH2O(J_0H:J_1H, LM, 12),
     *     ALB(IM, J_0H:J_1H, 9),
     *     SINJ(J_0H:J_1H),
     *     COSJ(J_0H:J_1H),
     &     srnflb_save(Im,J_0H:J_1H,Lm),
     &     trnflb_save(Im,J_0H:J_1H,Lm),
#ifdef TRACERS_ON
     &     ttausv_save(Im,J_0H:J_1H,Ntm,Lm),
     &     ttausv_cs_save(Im,J_0H:J_1H,Ntm,Lm),
#endif
     *     STAT=IER)

      KLIQ = 1
      dH2O = 0.
      SALB => ALB(:,:,1)
      SRVISSURF = 0
      RETURN
      END SUBROUTINE ALLOC_RAD_COM


      SUBROUTINE io_rad(kunit,iaction,ioerr) 3,57
!@sum  io_rad reads and writes radiation arrays to file
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : ioread,iowrite,irsfic,irerun,ioread_single
     *         ,lhead,Kradia,irsficnt,irsficno
#ifdef TRACERS_ON
      USE tracer_com,ONLY : Ntm
#endif
      USE RAD_COM
      USE PARAM
      USE DOMAIN_DECOMP, ONLY : GRID, GET, AM_I_ROOT
      USE DOMAIN_DECOMP, ONLY : UNPACK_COLUMN, PACK_COLUMN
      USE DOMAIN_DECOMP, ONLY : UNPACK_BLOCK , PACK_BLOCK
      USE DOMAIN_DECOMP, ONLY : UNPACK_DATA  , PACK_DATA
      USE DOMAIN_DECOMP, ONLY : ESMF_BCAST
      IMPLICIT NONE

      INTEGER kunit   !@var kunit unit number of read/write
      INTEGER iaction !@var iaction flag for reading or writing to file
!@var IOERR 1 (or -1) if there is (or is not) an error in i/o
      INTEGER, INTENT(INOUT) :: IOERR
!@var HEADER Character string label for individual records
      CHARACTER*80 :: HEADER, MODULE_HEADER = "RAD05"
!@var HEADER_F Character string label for records (forcing runs)
      CHARACTER*80 :: HEADER_F, MODULE_HEADER_F = "RADF"

      REAL*8 :: Tchg_glob(LM+LM_REQ,IM,JM)
      REAL*8 :: RQT_glob(LM_REQ,IM,JM)
      INTEGER :: kliq_glob(LM,4,IM,JM)
      REAL*8, DIMENSION(0:LM,IM,JM) :: SRHR_GLOB, TRHR_GLOB
      REAL*8, DIMENSION(1:4,IM,JM) :: TRSURF_GLOB
      REAL*8  :: FSF_GLOB(4,IM,JM)
      REAL*8, DIMENSION(IM, JM) :: FSRDIR_GLOB, SRVISSURF_GLOB,
     &           SRDN_GLOB, CFRAC_GLOB, SALB_GLOB
      REAL*8, DIMENSION(LM, IM, JM) :: RCLD_GLOB,O3_tracer_save_GLOB
      REAL*8, DIMENSION(LM, IM, JM, 5) :: rad_to_chem_GLOB
      REAL*8,DIMENSION(Im,Jm,Lm) :: srnflb_save_glob,trnflb_save_glob
#ifdef TRACERS_ON
      REAL*8,DIMENSION(Im,Jm,Ntm,Lm) :: ttausv_save_glob,
     &     ttausv_cs_save_glob
#endif
      INTEGER :: J_0,J_1
      INTEGER :: k

      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

      if (kradia.gt.0) then
        write (MODULE_HEADER_F(lhead+1:80),'(a8,i2,a15,i2,a7)')
     *    'R8 Tchg(',lm+LM_REQ,',ijm), I4 KLIQ(',lm,',4,ijm)'
        SELECT CASE (IACTION)
        CASE (:IOWRITE)            ! output to standard restart file
          CALL PACK_COLUMN(grid, Tchg, Tchg_glob)
          Call PACK_BLOCK( grid, kliq, kliq_glob)
          IF (AM_I_ROOT())
     &       WRITE (kunit,err=10) MODULE_HEADER_F,Tchg_glob,kliq_glob
        CASE (IOREAD:)
          SELECT CASE  (IACTION)
          CASE (ioread,irerun,ioread_single)  ! input for restart
            if (AM_I_ROOT() )
     &        READ (kunit,err=10) HEADER_F,Tchg_glob,kliq_glob
            CALL UNPACK_COLUMN(grid, Tchg_glob, Tchg)
            CALL UNPACK_BLOCK( grid, kliq_glob, kliq)
          CASE (IRSFIC)  ! only way to start frc. runs
            if (AM_I_ROOT() )
     &        READ (kunit,err=10) HEADER,RQT_glob,kliq_glob ; Tchg = 0.
            CALL UNPACK_COLUMN(grid, RQT_glob, RQT)
            CALL UNPACK_BLOCK( grid, kliq_glob, kliq)
            call sync_param( "Ikliq", Ikliq )
            if(Ikliq.eq.1) kliq=1  ! hysteresis init: equilibrium
            if(Ikliq.eq.0) kliq=0  ! hysteresis init: dry
          END SELECT
        END SELECT
      else

      MODULE_HEADER(lhead+1:80) = 'R8:Tq(3,ijm), I:LIQ(lm,4,ijm),'//
     *  ' S0,s+tHr,fs(4),fd,sv,al,sd,cf,o3'

      SELECT CASE (IACTION)
      CASE (:IOWRITE)            ! output to standard restart file
        CALL PACK_COLUMN(grid,  RQT,  RQT_glob)
        Call PACK_BLOCK( grid, kliq, kliq_glob)
        CALL PACK_COLUMN(grid, SRHR, SRHR_GLOB)
        CALL PACK_COLUMN(grid, TRHR, TRHR_GLOB)
        CALL PACK_COLUMN(grid, TRSURF, TRSURF_GLOB)
        CALL PACK_COLUMN(grid,  FSF,  FSF_GLOB)

        CALL PACK_DATA( grid, SALB     , SALB_GLOB)
        CALL PACK_DATA( grid, FSRDIR   , FSRDIR_GLOB)
        CALL PACK_DATA( grid, SRVISSURF, SRVISSURF_GLOB)
        CALL PACK_DATA( grid, SRDN     ,   SRDN_GLOB)
        CALL PACK_DATA( grid, CFRAC    ,  CFRAC_GLOB)

        CALL PACK_COLUMN(grid, RCLD          , RCLD_GLOB)
#ifdef TRACERS_SPECIAL_Shindell
        CALL PACK_COLUMN(grid, O3_tracer_save, O3_tracer_save_GLOB)
        CALL PACK_COLUMN(grid, rad_to_chem, rad_to_chem_GLOB)
#endif
        CALL PACK_DATA(grid,srnflb_save,srnflb_save_glob)
        CALL PACK_DATA(grid,trnflb_save,trnflb_save_glob)
#ifdef TRACERS_ON
        CALL PACK_DATA(grid,ttausv_save,ttausv_save_glob)
        CALL PACK_DATA(grid,ttausv_cs_save,ttausv_cs_save_glob)
#endif

        IF (AM_I_ROOT())
     *     WRITE (kunit,err=10) MODULE_HEADER,RQT_GLOB,KLIQ_GLOB
  ! rest needed only if MODRAD>0 at restart
     *      ,S0,  SRHR_GLOB, TRHR_GLOB, FSF_GLOB , FSRDIR_GLOB
     *       ,SRVISSURF_GLOB, SALB_GLOB, SRDN_GLOB, CFRAC_GLOB,
     *       RCLD_GLOB,  TRSURF_GLOB
#ifdef TRACERS_SPECIAL_Shindell
     *      ,O3_tracer_save_GLOB,rad_to_chem_GLOB
#endif
#ifdef TRACERS_DUST
     &      ,srnflb_save_glob,trnflb_save_glob,ttausv_save_glob
     &      ,ttausv_cs_save_glob
#endif
      CASE (IOREAD:)
        SELECT CASE  (IACTION)
        CASE (ioread,IRERUN)  ! input for restart, rerun or extension
          if (AM_I_ROOT() ) then
            READ (kunit,err=10) HEADER,RQT_glob,KLIQ_glob
  ! rest needed only if MODRAD>0 at restart
  ! rest needed only if MODRAD>0 at restart
     *       ,S0,  SRHR_GLOB, TRHR_GLOB, FSF_GLOB , FSRDIR_GLOB
     *       ,SRVISSURF_GLOB, SALB_GLOB, SRDN_GLOB, CFRAC_GLOB,
     *        RCLD_GLOB ,TRSURF_GLOB
#ifdef TRACERS_SPECIAL_Shindell
     *       ,O3_tracer_save_GLOB,rad_to_chem_GLOB
#endif
#ifdef TRACERS_DUST
     &       ,srnflb_save_glob,trnflb_save_glob,ttausv_save_glob
     &       ,ttausv_cs_save_glob
#endif
            IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
              PRINT*,"Discrepancy in module version ",HEADER,
     *               MODULE_HEADER
              GO TO 10
            END IF
          end if

          CALL ESMF_BCAST(grid, S0)
          CALL UNPACK_COLUMN(grid,  RQT_glob,  RQT)
          Call UNPACK_BLOCK( grid, kliq_glob, kliq)
          CALL UNPACK_COLUMN(grid, SRHR_glob, SRHR)
          CALL UNPACK_COLUMN(grid, TRHR_glob, TRHR)
          CALL UNPACK_COLUMN(grid, TRSURF_glob, TRSURF)
          CALL UNPACK_COLUMN(grid,  FSF_glob,  FSF)

          CALL UNPACK_DATA( grid, SALB_glob     , SALB)
          CALL UNPACK_DATA( grid, FSRDIR_glob   , FSRDIR)
          CALL UNPACK_DATA( grid, SRVISSURF_glob, SRVISSURF)
          CALL UNPACK_DATA( grid, SRDN_glob     ,   SRDN)
          CALL UNPACK_DATA( grid, CFRAC_glob    ,  CFRAC)

          CALL UNPACK_COLUMN(grid, RCLD_glob          ,RCLD)
#ifdef TRACERS_SPECIAL_Shindell
          CALL UNPACK_COLUMN(grid, O3_tracer_save_glob,
     &         O3_tracer_save)
          CALL UNPACK_COLUMN(grid, rad_to_chem_glob, rad_to_chem)
#endif
          CALL UNPACK_DATA(grid,srnflb_save_glob,srnflb_save)
          CALL UNPACK_DATA(grid,trnflb_save_glob,trnflb_save)
#ifdef TRACERS_ON
          CALL UNPACK_DATA(grid,ttausv_save_glob,ttausv_save)
          CALL UNPACK_DATA(grid,ttausv_cs_save_glob,ttausv_cs_save)
#endif


        CASE (IRSFIC,irsficnt,IRSFICNO)  ! restart file of prev. run
          if (AM_I_ROOT() )
     &      READ (kunit,err=10) HEADER,RQT_glob,KLIQ_glob
          CALL UNPACK_COLUMN(grid,  RQT_glob,  RQT)
          Call UNPACK_BLOCK( grid, kliq_glob, kliq)
          call sync_param( "Ikliq", Ikliq )
          if(Ikliq.eq.1) kliq=1  ! hysteresis init: equilibrium
          if(Ikliq.eq.0) kliq=0  ! hysteresis init: dry
        END SELECT
      END SELECT
      end if

      RETURN
 10   IOERR=1
      RETURN
      END SUBROUTINE io_rad