radaca.F90

Go to the documentation of this file.

!OPTIONS XOPT(HSFUN)
SUBROUTINE radaca &
 &( kidia , kfdia , klon , ktdia , klev &
 &, paprs , pgelam, psin  , pclon, pslon , pth &
 &, paer  , pozon &
 &)

!***********************************************************************
! CAUTION: THIS ROUTINE WORKS ONLY ON A NON-ROTATED, UNSTRETCHED GRID
!***********************************************************************

!**** *RADACA  - COMPUTES DISTRIBUTION OF AEROSOLS AND OZONE

!     PURPOSE.
!     --------

!**   INTERFACE.
!     ----------
!        CALL *RADACA* FROM *RADINT*

!        EXPLICIT ARGUMENTS :
!        --------------------
!     ==== INPUTS ===
!     ==== OUTPUTS ===

!        IMPLICIT ARGUMENTS :   NONE
!        --------------------

!     METHOD.
!     -------


!     EXTERNALS.
!     ----------

!          NONE

!     REFERENCE.
!     ----------

!        SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
!        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S"

!     AUTHOR.
!     -------
!     J.-J. MORCRETTE  E.C.M.W.F.    91/03/15

!     MODIFICATIONS.
!     --------------
!     J.-J. MORCRETTE  E.C.M.W.F.    93/03/15   OPERATIONAL CLIMATOLOGY
!     JJMorcrette 98-12-21 GISS volcanic aerosol climatology
!     JJMorcrette 99-09    monthly climatology of tropospheric aerosols

!     ADAPTATION TO zahir.idris.fr
!     ----------------------------
!     HGallée     04-11-17 Definition of the arguments 2 & 3 of LEGTRI
!-----------------------------------------------------------------------

#include "tsmbkind.h"

USE yomcst   , ONLY : r        ,rpi
USE yoeaerd  , ONLY : cvdaes   ,cvdael   ,cvdaeu   ,cvdaed   ,&
            &rcaeops  ,rcaeopl  ,rcaeopu  ,rcaeopd  ,rctrbga  ,&
            &rcvobga  ,rcstbga  ,rctrpt   ,raesc    ,raess    ,&
            &raelc    ,raels    ,raeuc    ,raeus    ,raedc    ,&
            &raeds
USE yoeozoc  , ONLY : cozqc    ,cozqs    ,cozhc    ,cozhs
USE yoerad   , ONLY : lhvolca  ,lnewaer
USE yoeaerc  , ONLY : rsinct   ,rsincv   ,repaer   ,&
            &rtaebc  ,rtaeor   ,rtaesd   ,rtaess   ,rtaesu   , &
            &rtaevo


IMPLICIT NONE


!     DUMMY INTEGER SCALARS
integer_m :: kfdia
integer_m :: kidia
integer_m :: klev
integer_m :: klon
integer_m :: ktdia

integer_m :: kcf, krint, kshift

integer_m :: kcp__radaca
integer_m :: kdim_radaca

!     -----------------------------------------------------------------

!*       0.1   ARGUMENTS.
!              ----------

real_b :: paprs(klon,klev+1), pgelam(klon), psin(klon)&
 &,   pclon(klon)       , pslon(klon) , pth(klon,klev+1)

real_b :: paer(klon,6,klev),pozon(klon,klev)
!     -----------------------------------------------------------------

!*       0.2   LOCAL ARRAYS.
!              -------------

integer_m :: iinla1(klon), iinla2(klon)
integer_m :: iinlo1(klon), iinlo2(klon)

real_b :: zaed(klon), zael(klon), zaes(klon), zaeu(klon)
real_b :: zaeqdn(klon), zaeqdo(klon), zaeqln(klon), zaeqlo(klon)
real_b :: zaeqsn(klon), zaeqso(klon), zaequn(klon), zaequo(klon)

real_b :: zaerbc(klon), zaeror(klon), zaersd(klon)
real_b :: zaerss(klon), zaersu(klon), zaervo(klon)

real_b :: zaetrn(klon),zaetro(klon)

real_b :: zalp(66)
real_b :: zdpn(klon)  , zdpo(klon)
real_b :: zfaed(21)    , zfael(21)    , zfaes(21)    , zfaeu(21)
real_b :: zfozq(11)    , zfozh(11)
real_b :: zgrth(klon)
real_b :: zlon(klon)  , zlonr(72)    , znlo1(klon) , znlo2(klon)

real_b :: zozh(klon) , zozq(klon)
real_b :: zqofn(klon) , zqofo(klon)
real_b :: zsilat(klon), zsinr(46)

!     LOCAL INTEGER SCALARS
integer_m :: il, imm, imnc, imns, inla, inla1, inla2, inlo1, inlo2, &
          &itotpt, jk, jl, jlr, jmm, &
          &jnn, nlatr, nlonr, jaer, jend, jil, jjl, iprint, itot

!     LOCAL REAL SCALARS
real_b :: zaetr, zcos1, zcos10, zcos2, zcos3, zcos4,&
          &zcos5, zcos6, zcos7, zcos8, zcos9, zcphn3, &
          &zcpho3, zdpnmo, zgridr, zlatr, zsdpn3, zsdpo3, &
          &zsin, zsin1, zsin10, zsin2, zsin3, zsin4, &
          &zsin5, zsin6, zsin7, zsin8, zsin9
real_b :: zaerbc1, zaerbc2, zaeror1, zaeror2, zaersd1, zaersd2, &
          &zaerss1, zaerss2, zaersu1, zaersu2
real_b :: zdlonr

!     ------------------------------------------------------------------

!*         1.     "NEW AEROSOL DISTRIBUTION" PARAMETERS COMPUTATIONS
!                 --------------------------------------------------


!*         1.1    VOLCANIC AEROSOL DISTRIBUTION PARAMETERS
!                 ----------------------------------------

!                 GISS CLIMATOLOGY 
!                 ----------------

kcf=0
kshift=0
krint=1

IF (lhvolca) THEN
  nlatr=46
  zgridr=180._jprb/(nlatr-1)
  DO jlr=1,nlatr
    zlatr=90._jprb-(jlr-1)*zgridr
    zsinr(jlr)=sin(zlatr*rpi/180._jprb)
  ENDDO

  il=kshift
  DO jl=kidia,kfdia
    il=il+1
    inla=0
    zsilat(il)=-9999._jprb
    zsin=psin(jl)
    DO jlr=1,nlatr-1
      IF (zsin <= zsinr(jlr) .AND. zsin > zsinr(jlr+1)) THEN
        inla=jlr
        zsilat(il)=(zsin-zsinr(inla))/(zsinr(inla+1)-zsinr(inla))
        zaervo(il)=rtaevo(inla)+zsilat(il)*(rtaevo(inla+1)-rtaevo(inla))
      ENDIF
    ENDDO
    IF (zsin <= zsinr(nlatr-1) .AND. zsin >= zsinr(nlatr))THEN
      inla=nlatr
      zsilat(il)=(zsin-zsinr(inla-1))/(zsinr(inla)-zsinr(inla-1))
      zaervo(il)=rtaevo(inla-1)&
       &+zsilat(il)*(rtaevo(inla)-rtaevo(inla-1))
    ENDIF
    IF (inla == 0) THEN
!      CALL ABOR1(' Problem with lat. interpolation in radaca!')
      stop ' Problem with lat. interpolation in radaca!'
    ENDIF
  ENDDO

!                 TANRE ET AL. CLIMATOLOGY 
!                 ------------------------

ELSE
  il = kshift
  DO jl=kidia,kfdia
    il = il+1
    zaervo(il)=rcvobga
  ENDDO
ENDIF
itotpt=il



!*         1.2    TROPOSPHERIC AEROSOL DISTRIBUTION PARAMETERS
!                 --------------------------------------------

IF (lnewaer) THEN
!  print *,'LNEWAER= ',LNEWAER
!  print *,' inputs        SINLAT            LONGITUDE'
!  DO JL=KIDIA,KFDIA
!    print 9001,JL,PSIN(JL),PGELAM(JL)
9001    format(1x,'RADACA ',1x,i5,1x,2e15.8)
!  END DO         

!-- latitude       
  nlatr=46
  zgridr=180._jprb/(nlatr-1)
  DO jlr=1,nlatr
    zlatr=90._jprb-(jlr-1)*zgridr
    zsinr(jlr)=sin(zlatr*rpi/180._jprb)
  END DO
  nlonr=72
  zdlonr=2._jprb*rpi/nlonr
  DO jlr=1,nlonr
    zlonr(jlr)=(jlr-1)*zdlonr
  END DO

!  print *,'Reference grid for Tegen climatology'        
!  print 9121,(ZSINR(JLR),JLR=1,NLATR)
9121     format(1x,'ZSINR ',8e15.7)
!  print 9122,(ZLONR(JLR),JLR=1,NLONR)
9122     format(1x,'ZLONR ',8e15.7)
           
  il=kshift
  DO jl=kidia,kfdia         
    il=il+1
    iinla1(il)=0
    iinla2(il)=0
    zsilat(il)=-9999._jprb
    zsin=psin(jl)
    DO jlr=1,nlatr-1
      IF (zsin <= zsinr(jlr) .AND. zsin > zsinr(jlr+1)) THEN
        inla=jlr
        iinla1(il)=jlr
        iinla2(il)=jlr+1
        zsilat(il)=(zsin-zsinr(inla))/(zsinr(inla+1)-zsinr(inla))
      ENDIF
    ENDDO
    IF (zsin <= zsinr(nlatr-1) .AND. zsin >= zsinr(nlatr))THEN
      inla=nlatr
      iinla1(il)=nlatr-1
      iinla2(il)=nlatr
      zsilat(il)=(zsin-zsinr(inla-1))/(zsinr(inla)-zsinr(inla-1))
    END IF    
    IF (inla.EQ.0) THEN 
!      CALL ABOR1(' Problem with lat. interpolation in radaca!')
      stop ' Problem with lat. interpolation in radaca!'
    ENDIF
!    print 9123,JL,IL,PSIN(JL),INLA,ZSINR(INLA),ZSILAT(IL)
9123     format(1x,'Interp.Latit.',2i4,f10.7,i4,2f10.7)          
  END DO  

!-- longitude
  il=kshift
  DO jl=kidia,kfdia
    il=il+1
    iinlo1(il)=0
    iinlo2(il)=0
    zlon(il)=-9999.
    DO jlr=1,71
      IF (pgelam(jl) < zlonr(jlr+1) .AND. pgelam(jl) >= zlonr(jlr)) &
     &  THEN
        iinlo1(il)=jlr
        iinlo2(il)=jlr+1
        znlo1(il)=zlonr(jlr)
        znlo2(il)=zlonr(jlr+1)
      END IF
    END DO      
    IF (pgelam(jl) >= zlonr(72)) THEN
      iinlo1(il)=72
      iinlo2(il)= 1
      znlo1(il)=zlonr(72)
!      ZNLO2(IL)=ZLONR(72)+2.*RPI
      znlo2(il)=zlonr(72)+zdlonr
    ENDIF
! Martin control
!    print 9124,JL,IL,PGELAM(JL),IINLO1(IL),IINLO2(IL),ZNLO1(IL),ZNLO2(IL)
!9124     format(1x,'Interp.Longit0.',2I4,F10.7,2I5,2F10.7)          
  END DO  
 
  il=kshift
  DO jl=kidia,kfdia        
    il=il+1
    IF (iinlo1(il).EQ.0 .OR. iinlo2(il).EQ.0) THEN 
!      CALL ABOR1(' Problem with long. interpolation in radaca!')
      stop ' Problem with long. interpolation in radaca!'
    ENDIF
    zlon(il)=(pgelam(jl)-znlo1(il))/(znlo2(il)-znlo1(il))
    inlo1=iinlo1(il)
    inlo2=iinlo2(il)
    inla1=iinla1(il)
    inla2=iinla2(il)
    
    zaerbc1=rtaebc(inlo1,inla1) &
     &      +zsilat(il)*(rtaebc(inlo1,inla2)-rtaebc(inlo1,inla1))  
    zaerbc2=rtaebc(inlo2,inla1) &
     &      +zsilat(il)*(rtaebc(inlo2,inla2)-rtaebc(inlo2,inla1))  
    zaerbc(il)=zaerbc1+zlon(il)*(zaerbc2-zaerbc1)
        
    zaeror1=rtaeor(inlo1,inla1) &
     &      +zsilat(il)*(rtaeor(inlo1,inla2)-rtaeor(inlo1,inla1))  
    zaeror2=rtaeor(inlo2,inla1) &
     &      +zsilat(il)*(rtaeor(inlo2,inla2)-rtaeor(inlo2,inla1))  
    zaeror(il)=zaeror1+zlon(il)*(zaeror2-zaeror1)
        
    zaersd1=rtaesd(inlo1,inla1) &
     &      +zsilat(il)*(rtaesd(inlo1,inla2)-rtaesd(inlo1,inla1))  
    zaersd2=rtaesd(inlo2,inla1) &
     &      +zsilat(il)*(rtaesd(inlo2,inla2)-rtaesd(inlo2,inla1))  
    zaersd(il)=zaersd1+zlon(il)*(zaersd2-zaersd1)
        
    zaerss1=rtaess(inlo1,inla1) &
     &      +zsilat(il)*(rtaess(inlo1,inla2)-rtaess(inlo1,inla1))  
    zaerss2=rtaess(inlo2,inla1) &
     &      +zsilat(il)*(rtaess(inlo2,inla2)-rtaess(inlo2,inla1))  
    zaerss(il)=zaerss1+zlon(il)*(zaerss2-zaerss1)
        
    zaersu1=rtaesu(inlo1,inla1) &
     &      +zsilat(il)*(rtaesu(inlo1,inla2)-rtaesu(inlo1,inla1))  
    zaersu2=rtaesu(inlo2,inla1) &
     &      +zsilat(il)*(rtaesu(inlo2,inla2)-rtaesu(inlo2,inla1))  
    zaersu(il)=zaersu1+zlon(il)*(zaersu2-zaersu1)
          
!    print 9125,JL,IL,PSIN(JL),PGELAM(JL),ZSILAT(IL) &
!     &      ,RTAESU(INLO1,INLA2),RTAESU(INLO1,INLA1),ZAERSU1 &
!     &      ,RTAESU(INLO2,INLA2),RTAESU(INLO2,INLA1),ZAERSU2 &
!     &                    ,INLA1,INLO1,INLO2,INLA2
!9125     format(1x,'Interp.Longit1.',2I4,9F10.7,4I5)          
!    print 9126,JL,IL,PSIN(JL),PGELAM(JL),ZSILAT(IL),ZLON(IL) &
!     &      ,ZNLO1(IL),ZNLO2(IL),INLA1,INLO1,INLO2,INLA2
!9126     format(1x,'Interp.Longit2.',2I4,6F10.7,4I5)          
!    print 9127,JL,IL,ZAERBC(IL),ZAEROR(IL),ZAERSD(IL),ZAERSS(IL) &
!     &      ,ZAERSU(IL)
!9127     format(1x,'Interp.Longit3.',2I4,5F10.7)          
  END DO       
END IF

!     ------------------------------------------------------------------   

!*       2.     OZONE
!               -----

zsin=psin(kidia)

!*       2.1     CALL TO LEGTRI.
!                ---------------
kcp__radaca= 6
kdim_radaca=66
!***
CALL legtri (zsin,kcp__radaca,kdim_radaca,zalp)
!***

!*       2.2     LEGENDRE TRANSFORM FOR OZONE.
!                -----------------------------

DO jmm=1,11
  zfozq(jmm)=_zero_
  zfozh(jmm)=_zero_
ENDDO
imm=0
imnc=0
imns=0
DO jmm=1,6
  imm=imm+1
  DO jnn=jmm,6
    imnc=imnc+1
    zfozq(imm)=zfozq(imm)+zalp(imnc)*cozqc(imnc)
    zfozh(imm)=zfozh(imm)+zalp(imnc)*cozhc(imnc)
  ENDDO
  IF(jmm /= 1) THEN
    imm=imm+1
    DO jnn=jmm,6
      imns=imns+1
      zfozq(imm)=zfozq(imm)+zalp(imns+6)*cozqs(imns)
      zfozh(imm)=zfozh(imm)+zalp(imns+6)*cozhs(imns)
    ENDDO
  ENDIF
ENDDO

!*       2.3     FOURIER TRANSFORM FOR OZONE.
!                ----------------------------

il=kshift
DO jl=kidia,kfdia
  il=il+1
  zcos1=pclon(jl)
  zsin1=pslon(jl)
  zcos2=zcos1*zcos1-zsin1*zsin1
  zsin2=zsin1*zcos1+zcos1*zsin1
  zcos3=zcos2*zcos1-zsin2*zsin1
  zsin3=zsin2*zcos1+zcos2*zsin1
  zcos4=zcos3*zcos1-zsin3*zsin1
  zsin4=zsin3*zcos1+zcos3*zsin1
  zcos5=zcos4*zcos1-zsin4*zsin1
  zsin5=zsin4*zcos1+zcos4*zsin1
  zozq(il)=&
   &zfozq(1)+_two_*(zfozq(2)*zcos1+zfozq(3)*zsin1+zfozq(4)*zcos2 &
   &+zfozq(5)*zsin2+zfozq(6)*zcos3+zfozq(7)*zsin3+zfozq(8)&
   &*zcos4+zfozq(9)*zsin4+zfozq(10)*zcos5+zfozq(11)*zsin5)
  zozh(il)=&
   &zfozh(1)+_two_*(zfozh(2)*zcos1+zfozh(3)*zsin1+zfozh(4)*zcos2 &
   &+zfozh(5)*zsin2+zfozh(6)*zcos3+zfozh(7)*zsin3+zfozh(8)&
   &*zcos4+zfozh(9)*zsin4+zfozh(10)*zcos5+zfozh(11)*zsin5)
  zozh(il)=sqrt(zozh(il))**3
ENDDO

!     ------------------------------------------------------------------

!       3.     AEROSOLS
!              --------
!***
!       3.1     CALL TO LEGTRI

kcp__radaca=11
kdim_radaca=66
!***
CALL legtri (zsin,kcp__radaca,kdim_radaca,zalp)
!***

!       3.2     LEGENDRE TRANSFORM FOR AEROSOLS
!               -------------------------------

DO jmm=1,21
  zfaes(jmm) = _zero_
  zfael(jmm) = _zero_
  zfaeu(jmm) = _zero_
  zfaed(jmm) = _zero_
ENDDO
imm  = 0
imnc = 0
imns = 0
DO jmm=1,11
  imm  = imm+1
  DO jnn=jmm,11
    imnc = imnc+1
    zfaes(imm) = zfaes(imm)+zalp(imnc)*raesc(imnc)
    zfael(imm) = zfael(imm)+zalp(imnc)*raelc(imnc)
    zfaeu(imm) = zfaeu(imm)+zalp(imnc)*raeuc(imnc)
    zfaed(imm) = zfaed(imm)+zalp(imnc)*raedc(imnc)
  ENDDO
  IF(jmm /= 1) THEN
    imm  = imm+1
    DO jnn=jmm,11
      imns = imns+1
      zfaes(imm) = zfaes(imm)+zalp(imns+11)*raess(imns)
      zfael(imm) = zfael(imm)+zalp(imns+11)*raels(imns)
      zfaeu(imm) = zfaeu(imm)+zalp(imns+11)*raeus(imns)
      zfaed(imm) = zfaed(imm)+zalp(imns+11)*raeds(imns)
    ENDDO
  ENDIF
ENDDO

!       3.3     FOURIER TRANSFORM FOR AEROSOLS
!               ------------------------------

il = kshift
DO jl=kidia,kfdia
  il = il+1
  zcos1    = pclon(jl)
  zsin1    = pslon(jl)
  zcos2    = zcos1*zcos1-zsin1*zsin1
  zsin2    = zsin1*zcos1+zcos1*zsin1
  zcos3    = zcos2*zcos1-zsin2*zsin1
  zsin3    = zsin2*zcos1+zcos2*zsin1
  zcos4    = zcos3*zcos1-zsin3*zsin1
  zsin4    = zsin3*zcos1+zcos3*zsin1
  zcos5    = zcos4*zcos1-zsin4*zsin1
  zsin5    = zsin4*zcos1+zcos4*zsin1
  zcos6    = zcos5*zcos1-zsin5*zsin1
  zsin6    = zsin5*zcos1+zcos5*zsin1
  zcos7    = zcos6*zcos1-zsin6*zsin1
  zsin7    = zsin6*zcos1+zcos6*zsin1
  zcos8    = zcos7*zcos1-zsin7*zsin1
  zsin8    = zsin7*zcos1+zcos7*zsin1
  zcos9    = zcos8*zcos1-zsin8*zsin1
  zsin9    = zsin8*zcos1+zcos8*zsin1
  zcos10   = zcos9*zcos1-zsin9*zsin1
  zsin10   = zsin9*zcos1+zcos9*zsin1
  zaes(il) = zfaes(1) + _two_*&
   &( zfaes(2)*zcos1  + zfaes(3)*zsin1  + zfaes(4)*zcos2 &
   &+ zfaes(5)*zsin2  + zfaes(6)*zcos3  + zfaes(7)*zsin3 &
   &+ zfaes(8)*zcos4  + zfaes(9)*zsin4  + zfaes(10)*zcos5 &
   &+ zfaes(11)*zsin5 + zfaes(12)*zcos6 + zfaes(13)*zsin6 &
   &+ zfaes(14)*zcos7 + zfaes(15)*zsin7 + zfaes(16)*zcos8 &
   &+ zfaes(17)*zsin8 + zfaes(18)*zcos9 + zfaes(19)*zsin9 &
   &+ zfaes(20)*zcos10+ zfaes(21)*zsin10                 )
  zael(il) = zfael(1) + _two_*&
   &( zfael(2)*zcos1  + zfael(3)*zsin1  + zfael(4)*zcos2 &
   &+ zfael(5)*zsin2  + zfael(6)*zcos3  + zfael(7)*zsin3 &
   &+ zfael(8)*zcos4  + zfael(9)*zsin4  + zfael(10)*zcos5 &
   &+ zfael(11)*zsin5 + zfael(12)*zcos6 + zfael(13)*zsin6 &
   &+ zfael(14)*zcos7 + zfael(15)*zsin7 + zfael(16)*zcos8 &
   &+ zfael(17)*zsin8 + zfael(18)*zcos9 + zfael(19)*zsin9 &
   &+ zfael(20)*zcos10+ zfael(21)*zsin10                 )
  zaeu(il) = zfaeu(1) + _two_*&
   &( zfaeu(2)*zcos1  + zfaeu(3)*zsin1  + zfaeu(4)*zcos2 &
   &+ zfaeu(5)*zsin2  + zfaeu(6)*zcos3  + zfaeu(7)*zsin3 &
   &+ zfaeu(8)*zcos4  + zfaeu(9)*zsin4  + zfaeu(10)*zcos5 &
   &+ zfaeu(11)*zsin5 + zfaeu(12)*zcos6 + zfaeu(13)*zsin6 &
   &+ zfaeu(14)*zcos7 + zfaeu(15)*zsin7 + zfaeu(16)*zcos8 &
   &+ zfaeu(17)*zsin8 + zfaeu(18)*zcos9 + zfaeu(19)*zsin9 &
   &+ zfaeu(20)*zcos10+ zfaeu(21)*zsin10                 )
  zaed(il) = zfaed(1) + _two_*&
   &( zfaed(2)*zcos1  + zfaed(3)*zsin1  + zfaed(4)*zcos2 &
   &+ zfaed(5)*zsin2  + zfaed(6)*zcos3  + zfaed(7)*zsin3 &
   &+ zfaed(8)*zcos4  + zfaed(9)*zsin4  + zfaed(10)*zcos5 &
   &+ zfaed(11)*zsin5 + zfaed(12)*zcos6 + zfaed(13)*zsin6 &
   &+ zfaed(14)*zcos7 + zfaed(15)*zsin7 + zfaed(16)*zcos8 &
   &+ zfaed(17)*zsin8 + zfaed(18)*zcos9 + zfaed(19)*zsin9 &
   &+ zfaed(20)*zcos10+ zfaed(21)*zsin10                 )
ENDDO


!     ------------------------------------------------------------------

!*       4.      VERTICAL DISTRIBUTION
!*               ---------------------


il=kshift
DO jl=kidia,kfdia
  il=il+1
  zdpo(il)=paprs(jl,1)
  zcpho3=paprs(jl,1)**3
  zsdpo3=sqrt(zcpho3)
  IF (lnewaer) THEN
    zaeqso(il)= zaerss(il)*cvdaes(1)
    zaeqlo(il)=(zaeror(il)+zaersu(il))*cvdael(1)
    zaequo(il)= zaerbc(il)*cvdaeu(1)
    zaeqdo(il)= zaersd(il)*cvdaed(1)
  ELSE
    zaeqso(il)=rcaeops*zaes(il)*cvdaes(1)
    zaeqlo(il)=rcaeopl*zael(il)*cvdael(1)
    zaequo(il)=rcaeopu*zaeu(il)*cvdaeu(1)
    zaeqdo(il)=rcaeopd*zaed(il)*cvdaed(1)
  END IF  
  zaetro(il)=_one_
  zqofo(il)=zozq(il)*zsdpo3 / (zsdpo3 + zozh(il))
ENDDO

DO jk=1,klev
  il=kshift
  IF (kcf == 0) THEN
    DO jl=kidia,kfdia
      il=il+1
      zgrth(il)= pth(jl,jk)/pth(jl,jk+1)
    ENDDO
  ELSEIF (kcf == 1) THEN
    DO jl=kidia,kfdia
      il=il+1
      zgrth(il)= pth(il,jk)/pth(il,jk+1)
    ENDDO
  ENDIF

  il=kshift
  DO jl=kidia,kfdia
    il=il+1
    zdpn(il)=paprs(jl,jk+1)
    zcphn3=paprs(jl,jk+1)**3
    zsdpn3=sqrt(zcphn3)
    IF (lnewaer) THEN
      zaeqsn(il)= zaerss(il)*cvdaes(jk+1)
      zaeqln(il)=(zaeror(il)+zaersu(il))*cvdael(jk+1)
      zaequn(il)= zaerbc(il)*cvdaeu(jk+1)
      zaeqdn(il)= zaersd(il)*cvdaed(jk+1)
    ELSE
      zaeqsn(il)=rcaeops*zaes(il)*cvdaes(jk+1)
      zaeqln(il)=rcaeopl*zael(il)*cvdael(jk+1)
      zaequn(il)=rcaeopu*zaeu(il)*cvdaeu(jk+1)
      zaeqdn(il)=rcaeopd*zaed(il)*cvdaed(jk+1)
    END IF  

    IF (_half_*(paprs(jl,jk)+paprs(jl,jk+1)) < 999._jprb) THEN
! for models with top above 10hPa
      zaetrn(il)=_one_
      zaetro(il)=_one_
    ELSE
      zaetrn(il)=zaetro(il)*(min(_one_,     zgrth(il)         ))**rctrpt
    ENDIF

    zaetr=sqrt(zaetrn(il)*zaetro(il))
    zqofn(il)=zozq(il)*zsdpn3/(zsdpn3+zozh(il))
    zdpnmo    =zdpn(il)-zdpo(il)
    paer(il,1,jk)=(_one_-zaetr)*(rctrbga*zdpnmo+ zaeqln(il)-zaeqlo(il))
    paer(il,2,jk)=(_one_-zaetr)*(zaeqsn(il)-zaeqso(il))
    paer(il,3,jk)=(_one_-zaetr)*(zaeqdn(il)-zaeqdo(il))
    paer(il,4,jk)=(_one_-zaetr)*(zaequn(il)-zaequo(il))
!old volc  PAER(IL,JK,5)= ZAETR * RCVOBGA*ZDPNMO
    paer(il,5,jk)=   zaetr  * zaervo(il) * zdpnmo
    paer(il,6,jk)=   zaetr  * rcstbga*zdpnmo
!old RH dependence          
!         AADS(IL,JK)=MAX(RCAEADM, (RCAEADK(1)*PAER(IL,1,JK)
!           + RCAEADK(2)*PAER(IL,2,JK)+RCAEADK(3)*PAER(IL,3,JK))/ZDPNMO)
    pozon(il,jk)=zqofn(il)-zqofo(il)
!**** **************************************************
!**** **************************************************
  ENDDO
  il=kshift
  DO jl=kidia,kfdia
    il=il+1
    zdpo(il)=zdpn(il)
    zqofo(il)=zqofn(il)

    zaeqso(il)=zaeqsn(il)
    zaeqlo(il)=zaeqln(il)
    zaequo(il)=zaequn(il)
    zaeqdo(il)=zaeqdn(il)
    zaetro(il)=zaetrn(il)
  ENDDO
  
  
!-- diagnostics in case of problem  
  DO jaer=1,6
    il=kshift
    DO jl=kidia,kfdia
      il=il+1
      paer(il,jaer,jk)=max(paer(il,jaer,jk),repaer)
    END DO
    itot=il
  END DO   
!--
   
ENDDO

!     ------------------------------------------------------------------

RETURN
END SUBROUTINE radaca