Skip to content
Snippets Groups Projects
radiations.f90 133 KiB
Newer Older
VIE Benoit's avatar
VIE Benoit committed
!MNH_LIC Copyright 1995-2022 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
!    ########################
     MODULE MODI_RADIATIONS   
!    ########################
!
CONTAINS
!
!   ############################################################################
    SUBROUTINE RADIATIONS (TPFILE,OCLEAR_SKY,OCLOUD_ONLY,                      &
               KCLEARCOL_TM1,HEFRADL,HEFRADI,HOPWSW,HOPISW,HOPWLW,HOPILW,      &
               PFUDG, KDLON, KFLEV, KRAD_DIAG, KFLUX, KRAD, KAER, KSWB_OLD,    &
               KSWB_MNH,KLWB_MNH, KSTATM,KRAD_COLNBR,PCOSZEN,PSEA, PCORSOL,    &
               PDIR_ALB, PSCA_ALB,PEMIS, PCLDFR, PCCO2, PTSRAD, PSTATM,        &
               PTHT, PRT, PPABST, POZON, PAER, PDST_WL, PAER_CLIM, PSVT,       &
               PDTHRAD, PSRFLWD, PSRFSWD_DIR,PSRFSWD_DIF, PRHODREF, PZZ,       &
               PRADEFF, PSWU, PSWD, PLWU,PLWD, PDTHRADSW, PDTHRADLW            )
!   ############################################################################
!
!!****  *RADIATIONS * - routine to call the SW and LW radiation calculations
!!
!!    PURPOSE
!!    -------
!!      The purpose of this routine is to prepare the temperature, water vapor
!!    liquid water, cloud fraction, ozone profiles for the ECMWF radiation
!!    calculations. There is a great number of available radiative fluxes in
!!    the output, but only the potential temperature radiative tendency and the
!!    SW and LW surface fluxes are provided in the output of the routine.
!!    Two simplified computations are available (switches OCLEAR_SKY and
!!    OCLOUD_ONLY). When OCLOUD_ONLY is .TRUE. the computations are performed
!!    for the cloudy columns only. Furthermore with OCLEAR_SKY being .TRUE.
!!    the clear sky columns are averaged and the computations are made for
!!    the cloudy columns plus a single ensemble-mean clear sky column.
!!
!!**  METHOD
!!    ------
!!      First the temperature, water vapor, liquid water, cloud fraction
!!    and  profile arrays are built using the current model fields and
!!    the standard atmosphere for the upper layer filling.
!!    The standard atmosphere is used between the levels IKUP and
!!    KFLEV where KFLEV is the number of vertical levels for the radiation 
!!    computations.    
!!    The aerosols optical thickness and the ozone fields come directly
!!    from ini_radiation step (climatlogies used) and are already defined for KFLEV. 
!!    Surface parameter ( albedo, emiss ) are also defined from current surface fields.
!!    In the case of clear-sky or cloud-only approximations, the cloudy
!!    columns are selected by testing the vertically integrated cloud fraction
!!    and the radiation computations are performed for these columns plus the
!!    mean clear-sky one. In addition, columns where cloud have disapeared are determined
!!    by saving cloud trace between radiation step and they are also recalculated
!!    in cloud only step. In all case, the sun position correponds to  the centered
!!    time between 2 full radiation steps (determined in physparam).
!!      Then the ECMWF radiation package is called and the radiative
!!    heating/cooling tendancies are reformatted in case of partial
!!    computations.  In case of "cloud-only approximation" the only cloudy
!!    column radiative fields are updated.
!!
!!    EXTERNAL
!!    --------
!!      Subroutine ECMWF_RADIATION_VERS2 : ECMWF interface calling radiation routines
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      Module MODD_CST  : constants
!!        XP00 : reference pressure
!!        XCPD : calorific capacity of dry air at constant pressure
!!        XRD  : gas constant for dry air
!!      Module MODD_PARAMETERS : parameters
!!        JPHEXT : Extra columns on the horizontal boundaries
!!        JPVEXT : Extra levels on the vertical boundaries
!!
!!    REFERENCE
!!    ---------
!!      Book2 of documentation ( routine RADIATIONS )
!!
!!    AUTHOR
!!    ------
!!	J.-P. Pinty      * Laboratoire d'Aerologie*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    26/02/95 
!!      J.Stein     20/12/95 add the array splitting in order to save memory
!!      J.-P. Pinty 19/11/96 change the split arrays, specific humidity
!!                           and add the ice phase
!!      J.Stein     22/06/97 use of the absolute pressure
!!      P.Jabouille 31/07/97 impose a zero humidity for dry simulation
!!      V.Masson    22/09/97 case of clear-sky approx. with no clear-sky column
!!      V.Masson    07/11/97 half level pressure defined from averaged Exner
!!                           function
!!      V.Masson    07/11/97 modification of junction between standard atm
!!                           and model for half level variables (top model
!!                           pressure and temperatures are used preferentially
!!                           to atm standard profile for the first point).
!!      P.Jabouille 24/08/98 impose positivity for ZQLAVE
!!      J.-P. Pinty 29/01/98 add storage for diagnostics
!!      J. Stein    18/07/99 add the ORAD_DIAG switch and keep inside the
!!                           subroutine the partial tendencies 
!!
!!      F.Solmon    04/03/01  MAJOR MODIFICATIONS, updated version of ECMWF radiation scheme
!!      P.Jabouille 05/05/03 bug in humidity conversion
!!      Y.Seity     25/08/03  KSWB=6 for SW direct and scattered surface 
!!                            downward fluxes used in surface scheme. 
!!      P. Tulet    01/20/05  climatologic SSA
!!      A. Grini    05/20/05  dust direct effect (optical properties)
!!      V.Masson, C.Lac 08/10 Correction of inversion of Diffuse and direct albedo
!!      B.Aouizerats 2010     Explicit aerosol optical properties
!!      C.Lac       11/2015   Correction on aerosols
!!      B.Vie            /13  LIMA
!!      J.Escobar 30/03/2017  : Management of compilation of ECMWF_RAD in REAL*8 with MNH_REAL=R4
!!      J.Escobar 29/06/2017  : Check if Pressure Decreasing with height <-> elsif PB & STOP 
!!      Q.LIBOIS  06/2017     : correction on CLOUD_ONLY
!!      Q.Libois  02/2018     : ECRAD
!!      Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!!      J.Escobar 28/06/2018 : Reproductible parallelisation of CLOUD_ONLY case
!!      J.Escobar 20/07/2018 : for real*4 compilation, convert with REAL(X) argument to SUM_DD... 
!!      P.Wautelet 22/01/2019: use standard FLUSH statement instead of non standard intrinsics
!!      Bielli S. 02/2019  Sea salt : significant sea wave height influences salt emission; 5 salt modes
!  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
!  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
!  P. Wautelet 06/09/2022: small fix: GSURF_CLOUD was not set outside of physical domain
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE PARKIND1,         ONLY: JPRB
USE OYOESW    , ONLY : RTAUA    ,RPIZA    ,RCGA
!
USE MODD_CH_AEROSOL,  ONLY: LORILAM
USE MODD_CONF,        ONLY: LCARTESIAN
USE MODD_CST
USE MODD_DUST,        ONLY: LDUST
use modd_field,          only: tfieldmetadata, TYPEREAL
VIE Benoit's avatar
VIE Benoit committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538
USE MODD_GRID ,       ONLY: XLAT0, XLON0
USE MODD_GRID_n ,     ONLY: XLAT, XLON
USE MODD_IO,          ONLY: TFILEDATA
USE MODD_LUNIT_n,     ONLY: TLUOUT
USE MODD_NSV,         ONLY: NSV_C2R2,NSV_C2R2BEG,NSV_C2R2END,     &
                            NSV_C1R3,NSV_C1R3BEG,NSV_C1R3END,     &
                            NSV_DSTBEG, NSV_DSTEND,               &
                            NSV_AERBEG, NSV_AEREND,               &
                            NSV_SLTBEG, NSV_SLTEND,               &
                            NSV_LIMA,NSV_LIMA_BEG,NSV_LIMA_END,   &
                            NSV_LIMA_NC, NSV_LIMA_NR, NSV_LIMA_NI
USE MODD_PARAMETERS
USE MODD_PARAM_LIMA
USE MODD_PARAM_n,     ONLY: CCLOUD, CRAD
USE MODD_PARAM_RAD_n, ONLY: CAOP
USE MODD_RAIN_ICE_DESCR
USE MODD_SALT,        ONLY: LSALT
USE MODD_TIME
!
USE MODE_DUSTOPT
USE MODE_IO_FIELD_WRITE, only: IO_Field_write
USE MODE_ll
use mode_msg
USE MODE_REPRO_SUM,      ONLY : SUM_DD_R2_R1_ll,SUM_DD_R1_ll
!
#ifdef MNH_PGI
USE MODE_PACK_PGI
#endif
USE MODE_SALTOPT
USE MODE_SUM_ll,          ONLY: MIN_ll
USE MODE_SUM2_ll,         ONLY: GMINLOC_ll
USE MODE_THERMO
!
USE MODI_AEROOPT_GET
USE MODI_ECMWF_RADIATION_VERS2
USE MODI_ECRAD_INTERFACE
USE MODD_VAR_ll,      ONLY: IP
!  
IMPLICIT NONE
!
!*       0.1   DECLARATIONS OF DUMMY ARGUMENTS :
!
TYPE(TFILEDATA),  INTENT(IN)         :: TPFILE    ! Output file
LOGICAL, INTENT(IN)                  :: OCLOUD_ONLY! flag for the cloud column
                                                   !    computations only
LOGICAL, INTENT(IN)                  :: OCLEAR_SKY ! 
INTEGER, INTENT(IN)                  :: KDLON   ! number of columns where the
                                                ! radiation calculations are
                                                !       performed
INTEGER, INTENT(IN)                  :: KFLEV   ! number of vertical levels
                                                !    where the radiation
                                                ! calculations are performed
INTEGER, INTENT(IN)                  :: KRAD_DIAG  ! index for the number of
                                                   !  fields in the output
INTEGER, INTENT(IN)                  :: KFLUX   ! number of top and ground 
                                                ! fluxes for the ZFLUX array
INTEGER, INTENT(IN)                  :: KRAD    ! number of satellite radiances
                                                ! for the ZRAD and ZRADCS arrays
INTEGER, INTENT(IN)                  :: KAER    ! number of AERosol classes

INTEGER, INTENT(IN)                  :: KSWB_OLD    ! number of SW band ECMWF 
INTEGER, INTENT(IN)                  :: KSWB_MNH    ! number of SW band ECRAD
INTEGER, INTENT(IN)                  :: KLWB_MNH    ! number of LW band ECRAD
INTEGER, INTENT(IN)                  :: KSTATM  ! index of the standard 
                                                ! atmosphere level just above
                                                !      the model top
INTEGER, INTENT(IN)                  :: KRAD_COLNBR ! factor by which the memory
                                                    ! is split
                                                    !
                                               !Choice of :             
CHARACTER (LEN=*), INTENT (IN)       :: HEFRADL ! 
CHARACTER (LEN=*), INTENT (IN)       :: HEFRADI ! 
CHARACTER (LEN=*), INTENT (IN)       :: HOPWSW !cloud water SW optical properties   
CHARACTER (LEN=*), INTENT (IN)       :: HOPISW !ice water SW optical properties 
CHARACTER (LEN=*), INTENT (IN)       :: HOPWLW !cloud water LW optical properties
CHARACTER (LEN=*), INTENT (IN)       :: HOPILW !ice water  LW optical properties
REAL,               INTENT(IN)       :: PFUDG  ! subgrid cloud inhomogenity factor
REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN ! COS(zenithal solar angle)
REAL,                     INTENT(IN) :: PCORSOL ! SOLar constant CORrection
REAL, DIMENSION(:,:),     INTENT(IN) :: PSEA    ! Land-sea mask
!
REAL, DIMENSION(:,:,:),   INTENT(IN) :: PDIR_ALB! Surface direct ALBedo
REAL, DIMENSION(:,:,:),   INTENT(IN) :: PSCA_ALB! Surface diffuse ALBedo
REAL, DIMENSION(:,:,:),   INTENT(IN) :: PEMIS   ! Surface IR EMISsivity
REAL, DIMENSION(:,:,:),   INTENT(IN) :: PCLDFR  ! CLouD FRaction
REAL,                     INTENT(IN) :: PCCO2   ! CO2 content
REAL, DIMENSION(:,:),     INTENT(IN) :: PTSRAD  ! RADiative Surface Temperature
REAL, DIMENSION(:,:),     INTENT(IN) :: PSTATM  ! selected standard atmosphere
!
REAL, DIMENSION(:,:,:),   INTENT(IN) :: PTHT    ! THeta at t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT     ! moist variables at t (humidity, cloud water, rain water, ice water)
REAL, DIMENSION(:,:,:),   INTENT(IN) :: PPABST  ! pressure at t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT    ! scalar variable ( C2R2 and C1R3  particle)
!
REAL, DIMENSION(:,:,:),   POINTER    :: POZON   ! OZONE field from clim.
REAL, DIMENSION(:,:,:,:), POINTER    :: PAER    ! AERosols optical thickness from clim. 
REAL, DIMENSION(:,:,:,:), POINTER    :: PDST_WL    ! AERosols Extinction by wavelength . 
REAL, DIMENSION(:,:,:,:), POINTER    :: PAER_CLIM    ! AERosols optical thickness from clim.
                                                ! note : the vertical dimension of 
                                                ! these fields include the "radiation levels"
                                                ! above domain top
                                                ! 
                                                 
REAL, DIMENSION(:,:,:), INTENT(IN)   :: PRHODREF ![kg/m3] air density
REAL, DIMENSION(:,:,:), INTENT(IN)   :: PZZ      ![m] height of layers

INTEGER, DIMENSION(:,:), INTENT(INOUT)  :: KCLEARCOL_TM1 ! trace of cloud/clear col
                                                         ! at the previous radiation step
!                                                 
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PDTHRAD ! THeta RADiative Tendancy
REAL, DIMENSION(:,:),     INTENT(INOUT) :: PSRFLWD ! Downward SuRFace LW Flux
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PSRFSWD_DIR ! Downward SuRFace SW Flux DIRect 
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PSRFSWD_DIF ! Downward SuRFace SW Flux DIFfuse 
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PSWU ! upward SW Flux 
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PSWD ! downward SW Flux 
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PLWU ! upward LW Flux 
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PLWD ! downward LW Flux 
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PDTHRADSW ! dthrad sw 
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PDTHRADLW !  dthradsw
REAL, DIMENSION(:,:,:),     INTENT(INOUT) :: PRADEFF ! effective radius
!
!
!*       0.2   DECLARATIONS OF LOCAL VARIABLES
!
LOGICAL                         :: GNOCL     ! .TRUE. when no cloud is present
                                             !     with OCLEAR_SKY .TRUE.
LOGICAL                         :: GAOP      ! .TRUE. when CAOP='EXPL'
LOGICAL, DIMENSION(KDLON,KFLEV) :: GCLOUD    ! .TRUE. for the cloudy columns
LOGICAL, DIMENSION(KFLEV,KDLON) :: GCLOUDT   ! transpose of the GCLOUD array
LOGICAL, DIMENSION(KDLON)       :: GCLEAR_2D ! .TRUE. for the clear-sky columns
LOGICAL, DIMENSION(KDLON,KFLEV) :: GCLEAR    ! .TRUE. for all the levels of the 
                                             !                clear-sky columns
LOGICAL, DIMENSION(KDLON,KSWB_MNH)  :: GCLEAR_SWB! .TRUE. for all the bands of the  
                                             !                clear-sky columns
INTEGER, DIMENSION(:), ALLOCATABLE :: ICLEAR_2D_TM1 !
!
INTEGER :: JI,JJ,JK,JK1,JK2,JKRAD,JALBS! loop indices
!
INTEGER :: IIB           ! I index value of the first inner mass point
INTEGER :: IJB           ! J index value of the first inner mass point
INTEGER :: IKB           ! K index value of the first inner mass point
INTEGER :: IIE           ! I index value of the last inner mass point
INTEGER :: IJE           ! J index value of the last inner mass point
INTEGER :: IKE           ! K index value of the last inner mass point
INTEGER :: IKU           ! array size for the third  index
INTEGER :: IIJ           ! reformatted array index
INTEGER :: IKSTAE        ! level number of the STAndard atmosphere array
INTEGER :: IKUP          ! vertical level above which STAndard atmosphere data
                         ! are filled in
!
INTEGER :: ICLEAR_COL    ! number of    clear-sky columns
INTEGER :: ICLOUD_COL    ! number of    cloudy    columns
INTEGER :: ICLOUD        ! number of levels corresponding of the cloudy columns
INTEGER :: IDIM          ! effective number of columns for which the radiation
                         ! code is run
INTEGER :: INIR          ! index corresponding to NIR fisrt band (in SW)
!
REAL, DIMENSION(:,:), ALLOCATABLE   :: ZTAVE    ! mean-layer temperature
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZTAVE_RAD    ! mean-layer temperature
REAL, DIMENSION(:,:), ALLOCATABLE   :: ZPAVE    ! mean-layer pressure
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZPAVE_RAD    ! mean-layer pressure
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQSAVE   ! saturation specific humidity
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQVAVE   ! mean-layer specific humidity
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQLAVE   ! Liquid water KG/KG
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQRAVE   ! Rain water  KG/KG
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQIAVE   ! Ice water Kg/KG
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQLWC   ! liquid water content kg/m3
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQRWC   ! Rain water  content kg/m3
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZQIWC   ! ice water content  kg/m3
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCFAVE   ! mean-layer cloud fraction
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZO3AVE   ! mean-layer ozone content 
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZPRES_HL ! half-level pressure
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZT_HL    ! half-level temperature
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZDPRES   ! layer pressure thickness
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCCT_C2R2! Cloud water Concentarion (C2R2)
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCRT_C2R2! Rain water Concentarion (C2R2)
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCIT_C1R3! Ice water Concentarion (C2R2)
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCCT_LIMA! Cloud water Concentration(LIMA)
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCRT_LIMA! Rain water Concentration(LIMA)
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCIT_LIMA! Ice water Concentration(LIMA)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE :: ZAER     ! aerosol optical thickness
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZALBP    ! spectral surface albedo for direct radiations
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZALBD    ! spectral surface albedo for diffuse radiations 
REAL(KIND=JPRB), DIMENSION (:,:),  ALLOCATABLE :: ZEMIS    ! surface LW  emissivity 
REAL(KIND=JPRB), DIMENSION (:,:), ALLOCATABLE  :: ZEMIW    ! surface LW  WINDOW emissivity
REAL(KIND=JPRB), DIMENSION(:), ALLOCATABLE     :: ZTS      ! reformatted surface PTSRAD array 
REAL(KIND=JPRB), DIMENSION(:), ALLOCATABLE     :: ZLSM     ! reformatted land sea mask
REAL(KIND=JPRB), DIMENSION(:),   ALLOCATABLE   :: ZRMU0    ! Reformatted ZMU0 array
REAL(KIND=JPRB)                     :: ZRII0    ! corrected solar constant
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDTLW    ! LW temperature tendency
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDTSW    ! SW temperature tendency
REAL, DIMENSION(:,:), ALLOCATABLE :: ZNFLW_CS ! CLEAR-SKY LW NET FLUXES
REAL, DIMENSION(:,:), ALLOCATABLE :: ZNFLW    ! TOTAL LW NET FLUXES
REAL, DIMENSION(:,:), ALLOCATABLE :: ZNFSW_CS ! CLEAR-SKY SW NET FLUXES
REAL, DIMENSION(:,:), ALLOCATABLE :: ZNFSW    ! TOTAL SW NET FLUXES
REAL, DIMENSION(:,:), ALLOCATABLE :: ZFLUX_TOP_GND_IRVISNIR ! Top and 
                                                            ! Ground radiative FLUXes
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_DOWN ! DowNward SW Flux profiles
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_UP   ! UPward   SW Flux profiles
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFLUX_LW      !          LW Flux profiles
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZDTLW_CS ! LW Clear-Sky temp. tendency
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZDTSW_CS ! SW Clear-Sky temp. tendency
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_TOP_GND_IRVISNIR_CS ! Top and
                                                  !  Ground Clear-Sky radiative FLUXes
REAL, DIMENSION(:,:), ALLOCATABLE   :: ZSFSWDIR !surface SW direct flux
REAL, DIMENSION(:,:), ALLOCATABLE   :: ZSFSWDIF !surface SW diffuse flux

REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_ALB_VIS, ZPLAN_ALB_NIR
                        ! PLANetary ALBedo in VISible, Near-InfraRed regions
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_TRA_VIS, ZPLAN_TRA_NIR
                        ! PLANetary TRANsmission in VISible, Near-InfraRed regions
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_ABS_VIS, ZPLAN_ABS_NIR
                        ! PLANetary ABSorption in VISible, Near-InfraRed regions
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZEFCL_LWD, ZEFCL_LWU
                        ! EFective  DOWNward and UPward LW nebulosity (equivalent emissivities)
                        ! undefined if RRTM is used for LW
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLWP, ZFIWP
                        ! Liquid and Ice Water Path
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZRADLP, ZRADIP
                        ! Cloud liquid water and ice effective radius
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZEFCL_RRTM, ZCLSW_TOTAL
                        ! effective LW nebulosity ( RRTM case) 
                        ! and SW CLoud fraction for mixed phase clouds
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTAU_TOTAL, ZOMEGA_TOTAL, ZCG_TOTAL
                        ! effective optical thickness, single scattering albedo
                        ! and asymetry factor for mixed phase clouds
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_DOWN_CS, ZFLUX_SW_UP_CS
                        ! Clear-Sky  DowNward and UPward   SW Flux profiles
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFLUX_LW_CS
                        ! Thicknes of the mesh
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZDZ
!
REAL, DIMENSION(KDLON,KFLEV) :: ZZDTSW ! SW diabatic heating
REAL, DIMENSION(KDLON,KFLEV) :: ZZDTLW ! LW diabatic heating
REAL, DIMENSION(KDLON)       :: ZZTGVIS! SW surface flux in the VIS band
REAL, DIMENSION(KDLON)       :: ZZTGNIR! SW surface flux in the NIR band
REAL, DIMENSION(KDLON)       :: ZZTGIR ! LW surface flux in the IR bands
REAL, DIMENSION(KDLON,SIZE(PSRFSWD_DIR,3)) :: ZZSFSWDIR
!                                      ! SW direct surface flux   
REAL, DIMENSION(KDLON,SIZE(PSRFSWD_DIR,3)) :: ZZSFSWDIF
!                                      ! SW diffuse surface flux   
!
REAL, DIMENSION(KDLON)       :: ZCLOUD ! vertically summed cloud fraction
!
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: ZEXNT ! Exner function
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2))    :: ZLWD    ! surface Downward LW flux
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PSRFSWD_DIR,3)) :: ZSWDDIR ! surface
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PSRFSWD_DIR,3)) :: ZSWDDIF ! surface Downward SW diffuse flux
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3),KSWB_OLD) :: ZPIZAZ ! Aerosols SSA
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3),KSWB_OLD) :: ZTAUAZ ! Aerosols Optical Detph
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3),KSWB_OLD) :: ZCGAZ  ! Aerosols Asymetric factor
REAL :: ZZTGVISC    ! downward surface SW flux (VIS band) for clear_sky
REAL :: ZZTGNIRC    ! downward surface SW flux (NIR band) for clear_sky
REAL :: ZZTGIRC     ! downward surface LW flux for clear_sky
REAL, DIMENSION(SIZE(PSRFSWD_DIR,3)) :: ZZSFSWDIRC
!                   ! downward surface SW direct flux for clear sky
REAL, DIMENSION(SIZE(PSRFSWD_DIR,3)) :: ZZSFSWDIFC
!                   ! downward surface SW diffuse flux for clear sky
REAL, DIMENSION(KFLEV) :: ZT_CLEAR  ! ensemble mean clear-sky temperature
REAL, DIMENSION(KFLEV) :: ZP_CLEAR  ! ensemble mean clear-sky temperature
REAL, DIMENSION(KFLEV) :: ZQV_CLEAR ! ensemble mean clear-sky specific humidity
REAL, DIMENSION(KFLEV) :: ZOZ_CLEAR ! ensemble mean clear-sky ozone
REAL, DIMENSION(KFLEV) :: ZHP_CLEAR ! ensemble mean clear-sky half-lev. pression
REAL, DIMENSION(KFLEV) :: ZHT_CLEAR ! ensemble mean clear-sky half-lev. temp.
REAL, DIMENSION(KFLEV) :: ZDP_CLEAR ! ensemble mean clear-sky pressure thickness
REAL, DIMENSION(KFLEV,KAER) :: ZAER_CLEAR  ! ensemble mean clear-sky aerosols optical thickness
REAL, DIMENSION(KSWB_MNH)       :: ZALBP_CLEAR ! ensemble mean clear-sky surface albedo (parallel)
REAL, DIMENSION(KSWB_MNH)       :: ZALBD_CLEAR ! ensemble mean clear-sky surface albedo (diffuse)
REAL                        :: ZEMIS_CLEAR ! ensemble mean clear-sky surface emissivity
REAL                        :: ZEMIW_CLEAR ! ensemble mean clear-sky LW window
REAL                        :: ZRMU0_CLEAR ! ensemble mean clear-sky MU0
REAL                        :: ZTS_CLEAR   ! ensemble mean clear-sky surface temperature.
REAL                        :: ZLSM_CLEAR  !  ensemble mean clear-sky land sea-mask  
REAL                        :: ZLAT_CLEAR,ZLON_CLEAR
!
!work arrays
REAL, DIMENSION(:),   ALLOCATABLE :: ZWORK1, ZWORK2, ZWORK3, ZWORK
REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK4, ZWORK1AER, ZWORK2AER, ZWORK_GRID
LOGICAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2)) :: ZWORKL
!
!  split arrays used to split the memory required by the ECMWF_radiation 
!  subroutine, the fields have the same meaning as their complete counterpart
!
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZALBP_SPLIT, ZALBD_SPLIT
REAL(KIND=JPRB), DIMENSION(:),     ALLOCATABLE :: ZEMIS_SPLIT, ZEMIW_SPLIT
REAL(KIND=JPRB), DIMENSION(:),     ALLOCATABLE :: ZRMU0_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZCFAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZO3AVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZT_HL_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZPRES_HL_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZTAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZPAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE :: ZAER_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZDPRES_SPLIT
REAL(KIND=JPRB), DIMENSION(:),     ALLOCATABLE :: ZLSM_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQVAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQSAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQLAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQIAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQRAVE_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQRWC_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQLWC_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZQIWC_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:),   ALLOCATABLE :: ZDZ_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCCT_C2R2_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCRT_C2R2_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCIT_C1R3_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCCT_LIMA_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCRT_LIMA_SPLIT
REAL(KIND=JPRB), DIMENSION(:,:), ALLOCATABLE   :: ZCIT_LIMA_SPLIT
REAL(KIND=JPRB), DIMENSION(:),     ALLOCATABLE :: ZTS_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZSFSWDIR_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZSFSWDIF_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZNFLW_CS_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZNFLW_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZNFSW_CS_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZNFSW_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZDTLW_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZDTSW_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_TOP_GND_IRVISNIR_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_DOWN_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_UP_SPLIT
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFLUX_LW_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZDTLW_CS_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZDTSW_CS_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_TOP_GND_IRVISNIR_CS_SPLIT
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_ALB_VIS_SPLIT
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_ALB_NIR_SPLIT
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_TRA_VIS_SPLIT
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_TRA_NIR_SPLIT
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_ABS_VIS_SPLIT
REAL, DIMENSION(:),     ALLOCATABLE :: ZPLAN_ABS_NIR_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZEFCL_LWD_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZEFCL_LWU_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLWP_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFIWP_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZRADLP_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZRADIP_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZEFCL_RRTM_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZCLSW_TOTAL_SPLIT
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTAU_TOTAL_SPLIT
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZOMEGA_TOTAL_SPLIT
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZCG_TOTAL_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_DOWN_CS_SPLIT
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZFLUX_SW_UP_CS_SPLIT
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFLUX_LW_CS_SPLIT
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZPIZA_EQ_TMP        !Single scattering albedo of aerosols (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZIR        !Real part of the aerosol refractive index(lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZII        !Imaginary part of the aerosol refractive index (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCGA_EQ_TMP         !Assymetry factor aerosols            (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZTAUREL_EQ_TMP      !tau/tau_{550} aerosols               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZPIZA_DST_TMP        !Single scattering albedo of dust (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCGA_DST_TMP         !Assymetry factor dust            (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZTAUREL_DST_TMP      !tau/tau_{550} dust               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZPIZA_AER_TMP        !Single scattering albedo of aerosol from ORILAM (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCGA_AER_TMP         !Assymetry factor aerosol from ORILAM            (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZTAUREL_AER_TMP      !tau/tau_{550} aerosol from ORILAM               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZPIZA_SLT_TMP        !Single scattering albedo of sea salt (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCGA_SLT_TMP         !Assymetry factor of sea salt            (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZTAUREL_SLT_TMP      !tau/tau_{550} of sea salt               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: PAER_AER      !tau/tau_{550} aerosol from ORILAM               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: PAER_SLT      !tau/tau_{550} sea salt               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: PAER_DST     !tau/tau_{550} dust               (lon,lat,lev,wvl)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTAU550_EQ_TMP      !tau/tau_{550} aerosols               (lon,lat,lev,wvl)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE   :: ZPIZA_EQ            !Single scattering albedo of aerosols (points,lev,wvl)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE   :: ZCGA_EQ             !Assymetry factor aerosols            (points,lev,wvl)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE   :: ZTAUREL_EQ          !tau/tau_{550} aerosols               (points,lev,wvl)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE   :: ZPIZA_EQ_SPLIT      !Single scattering albedo of aerosols (points,lev,wvl)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE   :: ZCGA_EQ_SPLIT       !Assymetry factor aerosols            (points,lev,wvl)
REAL(KIND=JPRB), DIMENSION(:,:,:), ALLOCATABLE   :: ZTAUREL_EQ_SPLIT    !tau/tau_{550} aerosols               (points,lev,wvl)
REAL, DIMENSION(KFLEV,KSWB_OLD)           :: ZPIZA_EQ_CLEAR      !Single scattering albedo of aerosols (lev,wvl)
REAL, DIMENSION(KFLEV,KSWB_OLD)           :: ZCGA_EQ_CLEAR       !Assymetry factor aerosols            (lev,wvl)
REAL, DIMENSION(KFLEV,KSWB_OLD)           :: ZTAUREL_EQ_CLEAR    !tau/tau_{550} aerosols               (lev,wvl)
INTEGER                               :: WVL_IDX              !Counter for wavelength

!
INTEGER  :: JI_SPLIT          ! loop on the split array
INTEGER  :: INUM_CALL         ! number of CALL of the radiation scheme
INTEGER  :: IDIM_EFF          ! effective number of air-columns to compute
INTEGER  :: IDIM_RESIDUE      ! number of remaining air-columns to compute
INTEGER  :: IBEG, IEND        ! auxiliary indices
!
!
REAL, DIMENSION(SIZE(PDTHRAD,1),SIZE(PDTHRAD,2),SIZE(PDTHRAD,3)) &
     :: ZDTRAD_LW! LW temperature tendency
REAL, DIMENSION(SIZE(PDTHRAD,1),SIZE(PDTHRAD,2),SIZE(PDTHRAD,3)) &
     :: ZDTRAD_SW! SW temperature tendency
INTEGER             :: ILUOUT       ! Logical unit number for output-listing
INTEGER             :: IRESP        ! Return code of FM routines
REAL, DIMENSION(SIZE(PDTHRAD,1),SIZE(PDTHRAD,2),SIZE(PDTHRAD,3)) &
     :: ZSTORE_3D, ZSTORE_3D2! 3D work array for storage
REAL, DIMENSION(SIZE(PDTHRAD,1),SIZE(PDTHRAD,2)) &
     :: ZSTORE_2D   ! 2D work array for storage!
INTEGER                         :: JBAND       ! Solar band index
CHARACTER (LEN=4), DIMENSION(KSWB_OLD) :: YBAND_NAME  ! Solar band name
CHARACTER (LEN=2)               :: YDIR        ! Type of the data field
!
INTEGER :: ISWB ! number of SW spectral bands (between radiations and surface schemes)
INTEGER :: JSWB ! loop on SW spectral bands
INTEGER :: JAE  ! loop on aerosol class
TYPE(TFIELDMeTaDATA) :: TZFIELD2D, TZFIELD3D
VIE Benoit's avatar
VIE Benoit committed
!
REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: ZDZPABST
REAL :: ZMINVAL
INTEGER, DIMENSION(3) :: IMINLOC
INTEGER :: IINFO_ll
LOGICAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2)) :: GCLOUD_SURF
!
REAL(KIND=JPRB), DIMENSION(:),   ALLOCATABLE :: ZLON,ZLAT
REAL(KIND=JPRB), DIMENSION(:),   ALLOCATABLE :: ZLON_SPLIT,ZLAT_SPLIT
!
INTEGER                            :: ICLEAR_COL_ll
INTEGER, DIMENSION(:), ALLOCATABLE :: INDEX_ICLEAR_COL
REAL, DIMENSION(KFLEV)             :: ZT_CLEAR_DD  ! ensemble mean clear-sky temperature
REAL                               :: ZCLEAR_COL_ll , ZDLON_ll
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!
!*       1.    COMPUTE DIMENSIONS OF ARRAYS AND OTHER INDICES
!              ----------------------------------------------
!
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)  ! this definition must be coherent with
                                      ! the one used in ini_radiations routine
IKU = SIZE(PTHT,3)
IKB = 1 + JPVEXT
IKE = IKU - JPVEXT
!
IKSTAE = SIZE(PSTATM,1)
IKUP   = IKE-JPVEXT+1
! 
ISWB   = SIZE(PSRFSWD_DIR,3)
!
!-------------------------------------------------------------------------------
!*       1.1   CHECK PRESSURE DECREASING
!              -------------------------
ZDZPABST(:,:,1:IKU-1) = PPABST(:,:,1:IKU-1) - PPABST(:,:,2:IKU)
ZDZPABST(:,:,IKU) = ZDZPABST(:,:,IKU-1)
!
ZMINVAL=MIN_ll(ZDZPABST,IINFO_ll)
!
IF ( ZMINVAL <= 0.0 ) THEN
   ILUOUT = TLUOUT%NLU
   IMINLOC=GMINLOC_ll( ZDZPABST )
   WRITE(ILUOUT,*) ' radiation.f90 STOP :: SOMETHING WRONG WITH PRESSURE , ZDZPABST <= 0.0 '  
   WRITE(ILUOUT,*) ' radiation :: ZDZPABST ', ZMINVAL,' located at ',   IMINLOC
   FLUSH(unit=ILUOUT)
   call Print_msg( NVERB_FATAL, 'GEN', 'RADIATIONS', 'something wrong with pressure: ZDZPABST <= 0.0' )

ENDIF
!------------------------------------------------------------------------------
ALLOCATE(ZLAT(KDLON))
ALLOCATE(ZLON(KDLON))
IF(LCARTESIAN) THEN
  ZLAT(:) = XLAT0*(XPI/180.)
  ZLON(:) = XLON0*(XPI/180.)
ELSE
  DO JJ=IJB,IJE
    DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZLAT(IIJ) =  XLAT(JI,JJ)*(XPI/180.)
        ZLON(IIJ) =  XLON(JI,JJ)*(XPI/180.)
    END DO
  END DO
END IF
!-------------------------------------------------------------------------------
!
!*       2.    INITIALIZES THE MEAN-LAYER VARIABLES
!              ------------------------------------
!
ZEXNT(:,:,:)= ( PPABST(:,:,:)/XP00 ) ** (XRD/XCPD)
!
! Columns where radiation is computed are put on a single line
ALLOCATE(ZTAVE(KDLON,KFLEV))
ALLOCATE(ZQVAVE(KDLON,KFLEV))
ALLOCATE(ZQLAVE(KDLON,KFLEV))
ALLOCATE(ZQIAVE(KDLON,KFLEV))
ALLOCATE(ZCFAVE(KDLON,KFLEV))
ALLOCATE(ZQRAVE(KDLON,KFLEV))
ALLOCATE(ZQLWC(KDLON,KFLEV))
ALLOCATE(ZQIWC(KDLON,KFLEV))
ALLOCATE(ZQRWC(KDLON,KFLEV))
ALLOCATE(ZDZ(KDLON,KFLEV))
!
ZQVAVE(:,:) = 0.0
ZQLAVE(:,:) = 0.0
ZQIAVE(:,:) = 0.0
ZQRAVE(:,:) = 0.0
ZCFAVE(:,:) = 0.0
ZQLWC(:,:) = 0.0
ZQIWC(:,:) = 0.0
ZQRWC(:,:) = 0.0
ZDZ(:,:)=0.0
!
!COMPUTE THE MESH SIZE
DO JK=IKB,IKE
  JKRAD = JK-JPVEXT
  DO JJ=IJB,IJE
    DO JI=IIB,IIE
      IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
      ZDZ(IIJ,JKRAD)  =  PZZ(JI,JJ,JK+1) - PZZ(JI,JJ,JK)
      ZTAVE(IIJ,JKRAD)  = PTHT(JI,JJ,JK)*ZEXNT(JI,JJ,JK) ! Conversion potential temperature -> actual temperature
    END DO
  END DO
END DO
!
!  Check if the humidity mixing ratio is available
!
IF( SIZE(PRT(:,:,:,:),4) >= 1 ) THEN
  DO JK=IKB,IKE
    JKRAD = JK-JPVEXT
    DO JJ=IJB,IJE
      DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZQVAVE(IIJ,JKRAD) =MAX(0., PRT(JI,JJ,JK,1))
      END DO
    END DO
  END DO
END IF
!
!  Check if the cloudwater mixing ratio is available
!
IF( SIZE(PRT(:,:,:,:),4) >= 2 ) THEN
  DO JK=IKB,IKE
    JKRAD = JK-JPVEXT
    DO JJ=IJB,IJE
      DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZQLAVE(IIJ,JKRAD) = MAX(0.,PRT(JI,JJ,JK,2))
        ZQLWC(IIJ,JKRAD) = MAX(0.,PRT(JI,JJ,JK,2)*PRHODREF(JI,JJ,JK))
        ZCFAVE(IIJ,JKRAD) = PCLDFR(JI,JJ,JK)
      END DO
    END DO
  END DO
END IF
!
!  Check if the rainwater mixing ratio is available
!
IF( SIZE(PRT(:,:,:,:),4) >= 3 ) THEN
  DO JK=IKB,IKE
    JKRAD = JK-JPVEXT
    DO JJ=IJB,IJE
      DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZQRWC(IIJ,JKRAD) = MAX(0.,PRT(JI,JJ,JK,3)*PRHODREF(JI,JJ,JK))
        ZQRAVE(IIJ,JKRAD) = MAX(0.,PRT(JI,JJ,JK,3))
      END DO
    END DO
  END DO
END IF
!
!  Check if the cloudice mixing ratio is available
!
IF( SIZE(PRT(:,:,:,:),4) >= 4 ) THEN
  DO JK=IKB,IKE
    JKRAD = JK-JPVEXT
    DO JJ=IJB,IJE
      DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZQIWC(IIJ,JKRAD) = MAX(0.,PRT(JI,JJ,JK,4)*PRHODREF(JI,JJ,JK))
!        ZQIAVE(IIJ,JKRAD) = MAX( PRT(JI,JJ,JK,4)-XRTMIN(4),0.0 )
        ZQIAVE(IIJ,JKRAD) = MAX( PRT(JI,JJ,JK,4),0.0 )
      END DO
    END DO
  END DO
END IF
!
!  Standard atmosphere extension
!
DO JK=IKUP,KFLEV
  JK1 = (KSTATM-1)+(JK-IKUP)
  JK2 = JK1+1
  ZTAVE(:,JK)  = 0.5*( PSTATM(JK1,3)+PSTATM(JK2,3) )
  ZQVAVE(:,JK) = 0.5*( PSTATM(JK1,5)/PSTATM(JK1,4)+   &
                 PSTATM(JK2,5)/PSTATM(JK2,4)    )
END DO
!
!        2.1 pronostic water concentation fields (C2R2 coupling) 
!
IF( NSV_C2R2 /= 0 ) THEN
  ALLOCATE (ZCCT_C2R2(KDLON, KFLEV))
  ALLOCATE (ZCRT_C2R2(KDLON, KFLEV))
  ZCCT_C2R2(:, :) = 0.
  ZCRT_C2R2 (:,:) = 0.
  DO JK=IKB,IKE
    JKRAD = JK-JPVEXT
    DO JJ=IJB,IJE
      DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZCCT_C2R2 (IIJ,JKRAD) = MAX(0.,PSVT(JI,JJ,JK,NSV_C2R2BEG+1))
        ZCRT_C2R2 (IIJ,JKRAD) = MAX(0.,PSVT(JI,JJ,JK,NSV_C2R2BEG+2))
      END DO
    END DO
  END DO
ELSE 
  ALLOCATE (ZCCT_C2R2(0,0))
  ALLOCATE (ZCRT_C2R2(0,0))
END IF
!
IF( NSV_C1R3 /= 0 ) THEN
  ALLOCATE (ZCIT_C1R3(KDLON, KFLEV))
  ZCIT_C1R3 (:,:) = 0.
  DO JK=IKB,IKE
    JKRAD = JK-JPVEXT
    DO JJ=IJB,IJE
      DO JI=IIB,IIE
        IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
        ZCIT_C1R3 (IIJ,JKRAD) = MAX(0.,PSVT(JI,JJ,JK,NSV_C1R3BEG))
      END DO
    END DO
  END DO
ELSE 
  ALLOCATE (ZCIT_C1R3(0,0))
END IF
!
!
!        2.1*bis pronostic water concentation fields (LIMA coupling) 
!
IF( CCLOUD == 'LIMA' ) THEN
   ALLOCATE (ZCCT_LIMA(KDLON, KFLEV))
   ALLOCATE (ZCRT_LIMA(KDLON, KFLEV))
   ALLOCATE (ZCIT_LIMA(KDLON, KFLEV))
   ZCCT_LIMA(:, :) = 0.
   ZCRT_LIMA (:,:) = 0.
   ZCIT_LIMA (:,:) = 0.
   DO JK=IKB,IKE
      JKRAD = JK-JPVEXT
      DO JJ=IJB,IJE
         DO JI=IIB,IIE
            IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
            IF (NMOM_C.GE.2) ZCCT_LIMA(IIJ,JKRAD) = MAX(0.,PSVT(JI,JJ,JK,NSV_LIMA_NC))
            IF (NMOM_R.GE.2) ZCRT_LIMA(IIJ,JKRAD) = MAX(0.,PSVT(JI,JJ,JK,NSV_LIMA_NR))
            IF (NMOM_I.GE.2) ZCIT_LIMA(IIJ,JKRAD) = MAX(0.,PSVT(JI,JJ,JK,NSV_LIMA_NI))
         END DO
      END DO
   END DO
END IF
!
!-------------------------------------------------------------------------------
!
!*       3.    INITIALIZES THE HALF-LEVEL VARIABLES
!  	           ------------------------------------
!
ALLOCATE(ZPRES_HL(KDLON,KFLEV+1))
ALLOCATE(ZT_HL(KDLON,KFLEV+1))
!
DO JK=IKB,IKE+1
  JKRAD = JK-JPVEXT
  DO JJ=IJB,IJE
    DO JI=IIB,IIE
      IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
      ZPRES_HL(IIJ,JKRAD) = XP00 * (0.5*(ZEXNT(JI,JJ,JK)+ZEXNT(JI,JJ,JK-1)))**(XCPD/XRD)
    END DO
  END DO
END DO

!  Standard atmosphere extension - pressure
!* begining at ikup+1 level allows to use a model domain higher than 50km
!
DO JK=IKUP+1,KFLEV+1
  JK1 = (KSTATM-1)+(JK-IKUP)
  ZPRES_HL(:,JK) = PSTATM(JK1,2)*100.0 ! mb -> Pa
END DO
!
!  Surface temperature at the first level
!  and surface radiative temperature
ALLOCATE(ZTS(KDLON))
!
DO JJ=IJB,IJE
  DO JI=IIB,IIE
    IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
    ZT_HL(IIJ,1) = PTSRAD(JI,JJ)
    ZTS(IIJ) = PTSRAD(JI,JJ)
  END DO
END DO
!
!  Temperature at half levels
!
ZT_HL(:,2:IKE-JPVEXT) = 0.5*(ZTAVE(:,1:IKE-JPVEXT-1)+ZTAVE(:,2:IKE-JPVEXT))
!
DO JJ=IJB,IJE
  DO JI=IIB,IIE
    IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
    ZT_HL(IIJ,IKE-JPVEXT+1)  =  0.5*PTHT(JI,JJ,IKE  )*ZEXNT(JI,JJ,IKE  ) &
         + 0.5*PTHT(JI,JJ,IKE+1)*ZEXNT(JI,JJ,IKE+1)
  END DO
END DO
!
!  Standard atmosphere extension - temperature
!* begining at ikup+1 level allows to use a model domain higher than 50km
!
DO JK=IKUP+1,KFLEV+1
  JK1 = (KSTATM-1)+(JK-IKUP)
  ZT_HL(:,JK) = PSTATM(JK1,3)
END DO
!
!mean layer pressure and layer differential pressure (from half level variables)
!
ALLOCATE(ZPAVE(KDLON,KFLEV))
ALLOCATE(ZDPRES(KDLON,KFLEV))
DO JKRAD=1,KFLEV
  ZPAVE(:,JKRAD)=0.5*(ZPRES_HL(:,JKRAD)+ZPRES_HL(:,JKRAD+1))
  ZDPRES(:,JKRAD)=ZPRES_HL(:,JKRAD)-ZPRES_HL(:,JKRAD+1)
END DO
!-----------------------------------------------------------------------
!*       4.    INITIALIZES THE AEROSOLS and OZONE PROFILES from climatology
!	           -------------------------------------------
!
!        4.1    AEROSOL optical thickness
! EXPL -> defined online, otherwise climatology
IF (CAOP=='EXPL') THEN
   GAOP = .TRUE.
ELSE
   GAOP = .FALSE.
ENDIF
!
IF (CAOP=='EXPL') THEN
   ALLOCATE(ZPIZA_EQ_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZCGA_EQ_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZTAUREL_EQ_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))

   ALLOCATE(ZPIZA_DST_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZCGA_DST_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZTAUREL_DST_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD)) 
   ALLOCATE(PAER_DST(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3)))

   ALLOCATE(ZPIZA_AER_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZCGA_AER_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZTAUREL_AER_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(PAER_AER(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3)))

   ALLOCATE(ZPIZA_SLT_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZCGA_SLT_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(ZTAUREL_SLT_TMP(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3),KSWB_OLD))
   ALLOCATE(PAER_SLT(SIZE(PAER,1),SIZE(PAER,2),SIZE(PAER,3)))
   

   ALLOCATE(ZII(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),KSWB_OLD))
   ALLOCATE(ZIR(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),KSWB_OLD))

  ZPIZA_EQ_TMP = 0.
  ZCGA_EQ_TMP = 0.
  ZTAUREL_EQ_TMP = 0.

  ZPIZA_DST_TMP = 0.
  ZCGA_DST_TMP = 0.
  ZTAUREL_DST_TMP = 0

  ZPIZA_SLT_TMP = 0.
  ZCGA_SLT_TMP = 0.
  ZTAUREL_SLT_TMP = 0

  ZPIZA_AER_TMP = 0.
  ZCGA_AER_TMP = 0.
  ZTAUREL_AER_TMP = 0

  PAER_DST=0.
  PAER_SLT=0.
  PAER_AER=0.
  
 IF (LORILAM) THEN
   CALL AEROOPT_GET(                             &
        PSVT(IIB:IIE,IJB:IJE,:,NSV_AERBEG:NSV_AEREND)        &  !I [ppv]  aerosols concentration
VIE Benoit's avatar
VIE Benoit committed
        ,PZZ(IIB:IIE,IJB:IJE,:)                   &  !I [m] height of layers
        ,PRHODREF(IIB:IIE,IJB:IJE,:)              &  !I [kg/m3] density of air
        ,ZPIZA_AER_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:)   &  !O [-] single scattering albedo of aerosols
        ,ZCGA_AER_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:)    &  !O [-] assymetry factor for aerosols
        ,ZTAUREL_AER_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:) &  !O [-] opt.depth(wvl=lambda)/opt.depth(wvl=550nm)
        ,PAER_AER(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT)            &  !O [-] optical depth of aerosols at wvl=550nm
        ,KSWB_OLD                                    &  !I |nbr] number of shortwave bands
        ,ZIR(IIB:IIE,IJB:IJE,:,:) &  !O [-] opt.depth(wvl=lambda)/opt.depth(wvl=550nm)
        ,ZII(IIB:IIE,IJB:IJE,:,:) &  !O [-] opt.depth(wvl=lambda)/opt.depth(wvl=550nm)
        )
 ENDIF
 IF(LDUST) THEN
   CALL DUSTOPT_GET(                             &
        PSVT(IIB:IIE,IJB:IJE,:,NSV_DSTBEG:NSV_DSTEND)        &  !I [ppv] Dust scalar concentration
VIE Benoit's avatar
VIE Benoit committed
        ,PZZ(IIB:IIE,IJB:IJE,:)                   &  !I [m] height of layers
        ,PRHODREF(IIB:IIE,IJB:IJE,:)              &  !I [kg/m3] density of air
        ,ZPIZA_DST_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:)   &  !O [-] single scattering albedo of dust
        ,ZCGA_DST_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:)    &  !O [-] assymetry factor for dust
        ,ZTAUREL_DST_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:) &  !O [-] opt.depth(wvl=lambda)/opt.depth(wvl=550nm)
        ,PAER_DST(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT)            &  !O [-] optical depth of dust at wvl=550nm
        ,KSWB_OLD                                   &  !I |nbr] number of shortwave bands
        )
   DO WVL_IDX=1,KSWB_OLD
     PDST_WL(:,:,:,WVL_IDX) = ZTAUREL_DST_TMP(:,:,:,WVL_IDX)* PAER(:,:,:,3)             
   ENDDO
 ENDIF
 IF(LSALT) THEN
   CALL SALTOPT_GET(                             &
        PSVT(IIB:IIE,IJB:IJE,:,NSV_SLTBEG:NSV_SLTEND)        &  !I [ppv] sea salt scalar concentration
VIE Benoit's avatar
VIE Benoit committed
        ,PZZ(IIB:IIE,IJB:IJE,:)                   &  !I [m] height of layers
        ,PRHODREF(IIB:IIE,IJB:IJE,:)              &  !I [kg/m3] density of air
        ,PTHT(IIB:IIE,IJB:IJE,:)                  &  !I [K] potential temperature
        ,PPABST(IIB:IIE,IJB:IJE,:)                &  !I [hPa] pressure
        ,PRT(IIB:IIE,IJB:IJE,:,:)                 &  !I [kg/kg] water mixing ratio
        ,ZPIZA_SLT_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:)   &  !O [-] single scattering albedo of sea salt
        ,ZCGA_SLT_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:)    &  !O [-] assymetry factor for sea salt
        ,ZTAUREL_SLT_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,:) &  !O [-] opt.depth(wvl=lambda)/opt.depth(wvl=550nm)
        ,PAER_SLT(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT)            &  !O [-] optical depth of sea salt at wvl=550nm
        ,KSWB_OLD                                    &  !I |nbr] number of shortwave bands
        )
 ENDIF

 ZTAUREL_EQ_TMP(:,:,:,:)=ZTAUREL_DST_TMP(:,:,:,:)+ZTAUREL_AER_TMP(:,:,:,:)+ZTAUREL_SLT_TMP(:,:,:,:)
 
 PAER(:,:,:,2)=PAER_SLT(:,:,:)
 PAER(:,:,:,3)=PAER_DST(:,:,:)
 PAER(:,:,:,4)=PAER_AER(:,:,:)


 WHERE (ZTAUREL_EQ_TMP(:,:,:,:).GT.0.0)
  ZPIZA_EQ_TMP(:,:,:,:)=(ZPIZA_DST_TMP(:,:,:,:)*ZTAUREL_DST_TMP(:,:,:,:)+&
                    ZPIZA_AER_TMP(:,:,:,:)*ZTAUREL_AER_TMP(:,:,:,:)+&
                    ZPIZA_SLT_TMP(:,:,:,:)*ZTAUREL_SLT_TMP(:,:,:,:))/&
                    ZTAUREL_EQ_TMP(:,:,:,:) 
 END WHERE
 WHERE ((ZTAUREL_EQ_TMP(:,:,:,:).GT.0.0).AND.(ZPIZA_EQ_TMP(:,:,:,:).GT.0.0))
  ZCGA_EQ_TMP(:,:,:,:)=(ZPIZA_DST_TMP(:,:,:,:)*ZTAUREL_DST_TMP(:,:,:,:)*ZCGA_DST_TMP(:,:,:,:)+&
                   ZPIZA_AER_TMP(:,:,:,:)*ZTAUREL_AER_TMP(:,:,:,:)*ZCGA_AER_TMP(:,:,:,:)+&
                   ZPIZA_SLT_TMP(:,:,:,:)*ZTAUREL_SLT_TMP(:,:,:,:)*ZCGA_SLT_TMP(:,:,:,:))/&
                   (ZTAUREL_EQ_TMP(:,:,:,:)*ZPIZA_EQ_TMP(:,:,:,:))
 END WHERE

 ZTAUREL_EQ_TMP(:,:,:,:)=max(1.E-8,ZTAUREL_EQ_TMP(:,:,:,:))
 ZCGA_EQ_TMP(:,:,:,:)=max(1.E-8,ZCGA_EQ_TMP(:,:,:,:))
 ZPIZA_EQ_TMP(:,:,:,:)=max(1.E-8,ZPIZA_EQ_TMP(:,:,:,:))
 PAER(:,:,:,3)=max(1.E-8,PAER(:,:,:,3))
 ZPIZA_EQ_TMP(:,:,:,:)=min(0.99,ZPIZA_EQ_TMP(:,:,:,:))


ENDIF      
!
! Computes SSA, optical depth and assymetry factor for clear sky (aerosols)
ZTAUAZ(:,:,:,:) = 0.
ZPIZAZ(:,:,:,:) = 0.
ZCGAZ(:,:,:,:)  = 0.
DO WVL_IDX=1,KSWB_OLD
 DO JAE=1,KAER
      !Special optical properties for dust
      IF (CAOP=='EXPL'.AND.(JAE==3)) THEN
      !Ponderation of aerosol optical in case of explicit optical factor
      !ti
        ZTAUAZ(IIB:IIE,IJB:IJE,IKB:IKE,WVL_IDX)= ZTAUAZ(IIB:IIE,IJB:IJE,IKB:IKE,WVL_IDX) + &
                                                 PAER(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,JAE) * &
                                       ZTAUREL_EQ_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,WVL_IDX) 
      !wi*ti
        ZPIZAZ(IIB:IIE,IJB:IJE,IKB:IKE,WVL_IDX)= ZPIZAZ(IIB:IIE,IJB:IJE,IKB:IKE,WVL_IDX) + &
                                  PAER(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,JAE)                * &
                                  ZTAUREL_EQ_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,WVL_IDX) * &
                                  ZPIZA_EQ_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,WVL_IDX)
      !wi*ti*gi
        ZCGAZ(IIB:IIE,IJB:IJE,IKB:IKE,WVL_IDX) = ZCGAZ(IIB:IIE,IJB:IJE,IKB:IKE,WVL_IDX) + &
                                 PAER(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,JAE)                * &
                                 ZTAUREL_EQ_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,WVL_IDX) * &
                                 ZPIZA_EQ_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,WVL_IDX)   * &
                                 ZCGA_EQ_TMP(IIB:IIE,IJB:IJE,IKB-JPVEXT:IKE-JPVEXT,WVL_IDX)
      ELSE

      !Ponderation of aerosol optical properties 
      !ti